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AN INVESTIGATION OF VOIDS FORMATION 

MECHANISMS AND THEIR EFFECTS ON FREEZE AND THAW 
PROCESS OF LITHIUM AND LITHIUM FLUORIDE 

ABSTRACT 

In the first part of the research, the mechanisms of void formation during the 
cooldown and freezing of lithium coolant within the primary loop of SP-100 type systems 
are investigated. These mechanisms are: (a) homogeneous nucleation, (b) heterogeneous 
nucleation, (c) normal segregation of helium gas dissolved in liquid lithium, and (d) 
shrinkage of lithium during freezing. To evaluate the void formation potential due to 
segregation, a numerical scheme that couples the freezing and mass diffusion processes 
in both the solid and liquid regions is developed. The results indicated that the formation 
of He bubbles is unlikely by either homogeneous or heterogeneous nucleation during the 
cooldown process. However, homogeneous nucleation of He bubbles following the 
segregation of dissolved He in liquid Lithium ahead of the solid-liquid interface is likely 
to occur. Results also show that total volume of He void is insignificant when compared 
to that of shrinkage voids. 

In viewing this, the subsequent research focuses on the effects of shrinkage void 
forming during freezing of lithium on subsequent thaw processes are investigated using a 
numerical scheme that is based on a single (solid/liquid) cell approach. The cases of 
lithium-fluoride are also investigated to show the effect of larger volume shrinkage upon 
freezing on the freeze and thaw processes. Results show that a void forming at the wall 
appreciably reduces the solid-liquid interface velocity, during both freeze and thaw, and 
causes a substantial rise in the wall temperature during thaw. However, in the case of Li, 
the maximum wall temperature was much lower than the melting temperature of PWC- 
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11, which is used as the structure material in the SP-100 system. Hence, it is concluded 
that a formation of hot spots is unlikely during the startup or restart of the SP-100 system. 
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1. INTRODUCTION 


1.1 Background 

The SP-100 Space Nuclear Power System, which is currently in the engineering 
development and testing phase, incorporates a fast flux, nuclear reactor which is cooled 
by a pumped liquid lithium (Li). The Li in the Primary Heat Transport System (PHTS), 
transports the reactor’s thermal power to the Power Conversion Assemblies (PCAs), 
where it is partially converted into electricity. The residual heat is transported by liquid 
Li in the secondary loops from the Thermoelectric-Electromagnetic (TEM) pumps and 
the PCAs to the system’s heat pipe radiator. 

The operating temperatures of the primary and secondary Li coolants at the 
reactor nominal thermal power of 2.4 MW, are about 1350 K to 1375 K, and 800 K to 
850 K, respectively. Because Li has a high melting point of -454 K, it will be frozen 
during launch and upon thawing and heating up to a nominal operating temperature of 
1350 K, Li will undergo about 20-25% increase in volume (see Figure 1). To minimize 
the stresses in the walls of the reactor containment and/or the piping system during 
startup in orbit, voids are introduced throughout the PHTS and each primary and 
secondary loop is equipped with a gas separator/accumulator and an accumulator with 
gas loaded bellows, respectively. 

Following a planned or an accidental shutdown, liquid Li in the primary and 
secondary coolant loops of the SP-100 system will eventually freeze. In the current SP- 
100 system design, the freezing process is controlled only by the natural law of cooling. 
The increase in the Li density during cooldown and freezing will stimulate formation of 
voids. Depending on the distribution of these voids, many different situations will arise 
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during the restart period of the SP-100 system. Homogeneously distributed voids within 
the solidified Li coolant, for example, would absorb the stress due to the volume increase 
of Li during the thaw process effectively. In this case, the stress in the walls of the 
coolant pipes, fuel pins and core container may not exceed those that cause deformation 
or failure of the materials. Dl posed voids, however, may lead to (a) mechanical failure 
of the coolant loop pipes and fuel elements due to excessive stresses, (b) thermal failure 
of structure or formation of hot spots during thaw. It is, therefore, important to 
understand their formation mechanisms and their effect on the system safety during 
rethaw. 

12 Void Formation Mechanisms 

In metal casting, volume shrinkage during freezing and bubble nucleation of 
dissolved gas atoms in the melt cause the formation of voids and results in a porous or 
leaky structure [Brick et al. (1977), Cahn & Haasen (1983)]. 

In the SP-100 system, helium (He) and tritium gas atoms are produced by 
neutron-lithium (n-Li) interactions in the reactor core. While tritium gas produced in the 
reactor core will readily diffuse through the niobium (Nb) walls to outer space, most of 
the He atoms will remain in the PHTS due to the very small diffusivity of He. McGhee, 
El-Genk and Rothrock (1989) have shown that Li coolant in the PHTS will be saturated 
with He gas within several days of operation and that the gas separators will only remove 
the amount of He gas in excess of the saturation concentration. It is, therefore, expected 
that He gas bubble might form, as in the metal casting process, during the cooldown and 
freeze processes of the Li coolant 
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Two well known mechanisms of gas bubble formation are: (1) homogeneous 
nucleation (HN), and (2) heterogeneous nucleation (HTN). HN refers to the appearance 
of a gas nucleus in the bulk of liquid, away from container walls or suspended particles, 
while HTN refers to the formation of a gas nucleus on surfaces or suspended impurities. 
The gas concentration required for HN of bubbles could be orders of magnitude higher 
than the saturation concentration in Li. On the other hand, in a system rich in suspended 
impurities, the gas concentration required for bubble nucleation by HTN in liquid Li 
could be as low as the saturation concentration. 

Another important mechanism of gas bubble formation during the freezing 
process is the segregation of dissolved gas atoms or molecules in a freezing medium. 
Segregation is a solute redistribution phenomenon during a freezing process. When the 
solubility limit of a solute in the solid phase is less than the solute concentration in the 
liquid phase, rejection of the solute by the S-L interface occurs as the freezing of the 
liquid progresses resulting in a buildup of the solute ahead of the S-L interface. Chalmers 
(1964) pointed out that in the case where gas is dissolved in a solidifying medium, 
segregation of gas in the liquid ahead of the solidification front can lead to the formation 
of gas voids. Figures 2a and 2b show gas void formation by segregation during the 
freezing processes of water and metal, respectively. In Chalmers (1964), seven different 
types of segregation are discussed including: (1) normal segregation, (2) grain boundary 
segregation, (3) cellular segregation, (4) dendritic segregation, (5) inverse segregation, 
(6) coring and intercrystalline segregation, and (7) gravity segregation. In microgravity, 
the last mechanism will not occur but the other six may occur, depending on the 
solidification process. The segregation of dissolved He atoms in liquid Li ahead of the 
solid-liquid (S-L) interface is expected to occur in the freezing process of Li coolant 
which may lead to the formation of gas bubbles during the process. The segregation 
mechanisms (2) - (6) require higher solidification velocities than that expected for normal 
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segregation, type (1). Because of the high heat capacity and latent heat of fusion of Li 
and the low heat rejection rate, since radiation is the only mode of heat rejection in the 
current design of the SP-100 system, the cooldown and freezing processes are expected to 
be very slow, hence normal segregation will be the more likely mechanism for formation 
of He gas bubbles in the SP-100 system. 

The formation of He gas bubbles by normal segregation will be accompanied by 
volume shrinkage of Li upon freezing stimulating the formation of shrinkage voids. 
These voids may occur either at the solid/wall interface or at the liquid/wall interface (see 
Figure 3). While the former will occur in the absence of adhesion between solidified Li 
and Nb wall, the latter is expected if liquid Li poorly wets Nb. 

The potential for He gas bubble formation by HN, HTN, and normal segregation 
is investigated and the formation of voids due to volume shrinkage is quantified in 
Section 2. In Section 3, the research focuses on quantifying the thermal effects of 
shrinkage voids during the freeze and thaw processes. 

1.3 Objectives 

This research investigates: (1) the void formation mechanisms and the potential 
for forming voids during the cooldown and freezing processes of Li coolant in the SP-100 
system after shutdown of the reactor in space, and (2) the thermal effect of void 
formation on the freeze and thaw processes of Li coolant Also, investigated are the 
thermal effects of void formation on both freezing and subsequent thawing of LiF, which 
is being considered for thermal energy storage in space application by NASA. Results of 
Li and LiF are compared and discussed. 
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(a) Bubbles in ice cube . 



(b) Gas bubbles in a metal. 


Figure 2. Gas Bubble Formation by Segregation during the 
Freezing Process of Water and Metal 
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(b) Shrinkage Void at Wall /Solid Interface 


Figure 3. Schematics of Void Formation in a Solidifying 
Lithium Coolant in a SP-100 System 
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2. VOID FORMATION MECHANISMS DURING THE COOLDOWN AND 
FREEZING PROCESSES OF LITHIUM COOLANT IN THE SP-100 SYSTEM 

2.1 Homogeneous and Heterogeneous Nudeatlon of He Gas Bubbles 
During the Cooldown Process of LI Coolant 

2.1.1 Theoretical Background 

In order for a gas nucleus to be thermodynamically stable, its radius should be 
larger than a critical value. Van Stralen (1979) has shown that the supersaturation, 
SShN. necessary for HN of a gas nucleus of critical radius can be represented by: 

SShn = ACyCga ■ to/r^L, (1) 

where, AC h = Cjj - C^. On die other hand, the supersaturation required for HTN of gas 
bubbles, SSHTN. may be expressed as (Cole 1974): 

SS HTN = SShn F(Geo, 6), (2) 

where F is a complex function of wetting angle, 6, and surface geometry. For a flat 
surface F is only a function of the liquid-solid wetting angle, 6 

m - 1 ' 3 . (3) 

Since, F varies between 0 and 1, SShtN can be as small as 0 (when liquid Li does not 
wet the container’s wall or when the liquid is rich in solid or gas impurities) and as large 
as SShn (when liquid Li perfectly wets the wall or the liquid is pure). An expression for 

r cr in equation (1) is (Cole 1974): 
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In this equation Ej can be reasonably taken to lie between 0 and the energy of formation 

of a gas cluster of radius r cr (0 < E(j < 4K&J/3). and J can be assumed to lie somewhere 
between 10^ and 10^ nuc./ra^-s for what is normally thought of as a gas ebullition 
process. It can be shown that r cr is approximately 1 nm, and is not a strong function of 
either J or E(j. Therefore, in the following assessment of the potential of He bubble 
nucleation during the cooldown process of liquid Li in SP-100 system, a value of r cr = 

10'^ m is used. The values of Csat are determined as functions of liquid Li pressure and 
temperature using the following correlation suggested by Slotnick et al. (1965): 

Csat = PL ex P('^*^ * 23.0) ^ 

The reduction of Li coolant volume during the cooldown process will be 
compensated for by volume increase in the gas separator/accumulator (GSA). If we 

assume an ideal gas system in the GSA, Pl in Equation (1) can be calculated from the 
ideal gas law: 

P L V _ P o v 0 

T " T e . (6) 

The bubble nucleation potential (BNP) by either HN or HTN is defined as the 
ratio of actual supersaturation of the dissolved gas to supersaturation required for HN or 
HTN: 

BNP = SS act /SSHN (or SSHTN). (7) 

He gas bubbles will nucleate when BNP is greater than or equal to unity. 

2. 1 .2 Results and Discussion 


In order to assess the potential of He gas bubble nucleation in the SP-100 PHTS, 
the initial He concentration in Li coolant, Co, is taken equal to saturation value, Csat. at 



the nominal operating condition in the SP-100 prior to reactor shutdown (see Table 1). 
Using the system conditions listed in Table 1, the variations of Pl, Cs and BNP with Li 
temperature are calculated and the results are shown in Figures 4 and 5. In order for 
either HN or HTN of He bubble to occur the BNP should be higher than or equal to unity. 
However, as Figure 5 demonstrates, the BNP values for HN are much smaller than unity, 
hence HN of He gas bubbles is not possible throughout the entire cooldown process of Li 
coolant in SP-100. Furthermore, it can be deduced from Equation (2) and Figure 5 that in 
order for He gas bubbles to form by HTN, the value of F(e) in Equation (3) needs to be as 
low as 2.5 x 10* 3 , requiring a wetting angle e = 150" (see Figure 6). Such high wetting 
angle is not generally expected in a liquid metal-metal substrate system as in the SP-100. 
Liquid Li is known to perfectly wet metallic surfaces. Also, dewetting does not occur 
upon cooling of U below the wetting temperature (673.5 K) once it is heated above the 
wetting temperature and maintained at the temperature for a while [Hoffman (1990)]. 
Therefore, it could be concluded that the formation of He bubbles by HTN during the 
cooldown process in the SP-100 PHTS is highly unlikely, provided all the system 
conditions are remaining as fabricated. 

2.2 Potential Formation of Intermetallic Compounds During the Cooldown 

Process of Li Coolant and Their Relaton to the Formation of He Gas 

Bubbles 

It has been concluded in the previous section that the formation of He gas bubbles 
during the cooldown process of Li coolant is highly unlikely, provided all the system 
conditions are remaining as fabricated. However, should the characteristics of the metal 
substrate change due to the formation of intermetallic compounds, such as metal hydride 
or metal nitride of thickness as little as a few atom radii, the likelihood of dewetting of Li 
and the formation of He bubbles at the wall by HTN increases. 
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2.2. 1 Potential Formation of Intermeta llic Compounds 


Because of their high compatibility with Li, Nb and Nb-l%Zr are used in the SP- 
100 system to fabricate the PHTS pipes and fuel cladding, respectively. As mentioned 
previously, the (n-Li) reactions in the reactor core produce tritium [McGhee, El-Genk and 
Rothrock (1989)]. Although the tritium’s high diffusivity enables it to readily diffuse out 
from the PHTS, a small quantities of tritium may remain in the PHTS prior to reactor 
shutdown. The diffusivity of tritium in Nb rapidly decreases as the temperature of Nb is 
reduced [Alefeld and Volki (1978)], hence increasing the likelihood of Nb-T reaction 
during the cooldown process of Li in the SP-100 PHTS and the formation of Nb-T 
compounds at the inner surface of the pipes wall. 

Another source for forming intermetallic compounds in SP-100 PHTS is nitrogen 

(N2) gas in the nuclear fuel elements. Because, N2 gas will be produced due to the 
dissociation of UN fuel during the operation of the SP-100 system. The amount of N2 
gas dissociated depends on the fuel temperature and bumup. Should the fuel cladding 
breach, N2 gas will easily leak from the fuel into the Li coolant. It is also expected that 
N2 gas will eventually diffuse through fuel cladding into Li coolant 

Table 2 shows that tritium and nitrogen readily react with Nb and form hydride 
and nitride at the reaction temperatures indicated. Since, the nominal operating 
temperature in the PHTS of the SP-100 system is about 1300 to 1350 K, formation of 
stable NbN is likely on the surface of fuel cladding or on the inner surfaces of the PHTS 
pipes. As the Li cools down after a reactor shutdown, within a temperature range of 673 
<T<1273 K, the formation of a mixture of NbN and Nb2N is expected. When the Li 
coolant temperature is below 973 K, the formation of LiT occurs and deposits saline LiT 
in the Li coolant. During the last stage of the cooldown process (T<573 K), the 
formation of NbT is expected. 
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Figure 4. Variation of System Pressure and Saturation 

Concentration of He in Li with Lithium Temperature 
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Saturation Concentration of He in Li ( 10 ' 7 kg/m 





Nucleation Potential (BNP) 






F(Rs/^vn'^) 


HETEROGENEOUS NUCLEATION 
GEOMETRIC FACTOR 


SUPERSATURATION FOR HETEROGENEOUS 
NUCLEATION EQUALS F FACTOR TIMES 
SUPERSATURATION REQUIRED FOR 
FOR HOMOGENEOUS NUCLEATION. 



WETTING ANGLE 0 (Degrees) 


Figure 6. Heterogeneous Nucleation Geometric Factor 
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2.2.2 Discussion and Suggestion for FutureJnyestigation 


The deposit of LIT and the formation of NbT, NbN, and Nb2N on inner surface of 
PHTS walls will increase the possibility of dewetting of Li from the walls, hence 
increasing the potential of void formation by HTN. However, in order to evaluate the 
potential, wetting angles of Li and mixture of Li and LiT on NbT, NbN, and Nb2N are 
needed. Our literature search indicated that such information is not available, hence the 
potential of He gas bubble formation by HTN in conjunction with the formation and 
deposit of intermetallic compounds on the inner walls of coolant pipes, reactor core 
container, and the outer surface of fuel pins remain inconclusive. Experimental efforts to 
determine the wetting angles of Li on Nb-nitride and Nb-hydride substrates are, therefore, 
highly recommended. 

2.3 He Gas Bubble Formation by Segregation of He Gas Atoms During the 

Freezing Process of Li Coolant 

In the normal segregation process, the motion of He gas (solute) out of the 
solidified Li is in the direction of movement of solidification interface. The rate of He 
buildup ahead of the S-L interface is determined by the rate of He rejection from the solid 
at the solidification front and the rate of He diffusion away from the S-L interface into the 
liquid. The former process will tend to raise He concentration at the S-L interface, while 
the latter will reduce it The rate of He rejection, J(8(t)), can be expressed as: 

J(8(t» « (C(8, t) - C* ) V(t) . (8) 

The solidification velocity, V(t)( = d8(t)/dt), is obtained from the thermal analysis of the 
freezing process. In order to determine the potential of He gas formation by normal 
segregation, a physical model is developed in which the freezing of Li and the mass 
diffusion of He away from the S-L interface are coupled. In this model the void 
formation due to shrinkage of Li during freezing is also simulated. 
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2.3.1 Physical Model. 

For simplicity, a one dimensional system initially filled with liquid Li at its 
melting temperature is considered. The liquid Li is cooled by thermal radiation into 
space at one wall, x = 0, and is thermally insulated at the other wall, x = L. In Figure 3b, 
the shrinkage void forms at the wall/solidified crust interface, while in Figure 3a, it forms 
in the liquid at the insulated wall. In this chapter the model in Figure 3a is analyzed. The 
analysis of the model in Figure 3b will be described in chapter 3 in conjunction with the 
effect of shrinkage void formation on the freeze and thaw processes. In Figure 3a, the 
location and velocity of the boundary “B” of the liquid phase can be represented by Xr = 


L - rs(t), v(t) = - r(d5(t)/dt), respectively, where r = (ps/pl - 1). The He concentration in 
liquid Li is assumed uniform prior to the initiation of the freezing process and equal to its 
saturation concentration in liquid Li at 1350 K (Co = 4.26 x 10-6 kg/m3). The 
segregation model consists of two coupled models: (a) freezing model of liquid Li, and 
(b) mass diffusion model of He gas in liquid Li ahead of the S-L interface. Because the 
problem in Figure 3a entails two moving boundaries, S-L interface “A” and shrinkage 
void/liquid interface “B”, and a boundary condition of the third kind (radiation), it is 
highly nonlinear. Hence, an exact analytical solution is not attainable. Instead, a 
numerical approach is used to solve coupled freezing and gas diffusion processes. 
Freezing Model of Liquid Li 

For the one dimensional system shown in Figures 3a and 7a, the heat balance 
equations in the liquid and solid Li and initial and boundary conditions can be written as: 
(A) Liq ui d Li ( S^^Xr) 


(B) 




Solid Li (0 <; x <; S) 


3T, 




dxpdx 


(9) 
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( 10 ) 


(11a) 


(C) Initial and Boundary Conditions 
T(x,0) = T f , 


k^O, t) = ec{t 4 (0, t) - Tt)= ^(mt) - T a ) 


(lib) 


fj<*,t) =0 

3x 

T s (5,t) = T 1 j(5,t) = T ff 


(lie) 

(lid) 



(lie) 


Because of the high thermal conductivity of Li, the temperature differential in the 
frozen crust is expected to be small. Hence, constant properties can be assumed. When 
constant properties in each spatial mesh, Ax, within the solid and liquid Li are assumed, 
Equation 9 can be written in terms of enthalpy as [Shamsunder and Sparrow (1975), 
Hunter and Kuttler (1989)): 


af x+A * 

3tJ* 


pHdx + 


pvH - k— 
dx 


x+Ax 


= 0 


J x 


( 12 ) 


The dimensionless nodal enthalpy and temperature are defined as: 


0E 


-i- p 2LM*. 

PAxJax X 


a , E sE iIl 

X 


(13) 


It has been shown by Shamsunder and Sparrow (1975), that 6 and <|> have the following 
relationships: 

= a + bG, 

where, 

a = 0, b = 1 .for 0 < 0 (14) 

a = 0, b = 0 , for OS 0S 1. 
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Figure 7. Schematic Representation of Freezing 
and He Segregation Processes 
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Using these dimensionless parameters. Equation (9) is transformed into: 


3£ = j 

jj^jx+Ax y ae 


dt A: 

Kp Cp dx x 3x 

(15) 


It has been shown by Shamsunder and Sparrow (1975) that in the absence of 
convection term due to volume shrinkage. Equation (15) satisfies Equation (lie). It can 
also be shown that the same approach can be applied in the case with convection due to 
volume shrinkage to show that the energy conservation at the S-L interface is satisfied by 
Equation (15). Using backward finite difference formulation for the time derivative and 
forward difference formulation of the spacial derivatives Equation (15) is discretized as: 

y-r‘ - ■ [/■a tei'-tf-) . /jl) ( tr-«a> L 

At AxftL|C p J. Ax VCpjw Ax J At Ax (16) 

The heat balance equation for the solid [Equation (10)] can be transformed and 
linearized as in Equations (15) and (16) by noting that v in the solid region is zero in and 
the last term of the Equation (16) drops out The location of S-L interface, 6, is 

determined by 6 = xj - 6jAxj when the value of 0 of the i-th node is less than 1 and greater 
than 0, since the frozen fraction of the node, a, can be conveniently related to 0 as a = 1 - 
0 from the definition of 0 shown in Equation (13). The transformation and linearization 
of the initial and boundary conditions in Equation (11) are carried out consistently with 
those of the governing equations and the values of 0 and 4> at each time step are 
determined by a fully implicit scheme suggested by Shamsunder and Sparrow (1975). 
The program developed in accordance with the foregoing considerations 
(COPHASE.FOR) is listed in Appendix A-l and input data required to run the program 
(COPHASE.INP) is listed in Appendix A-2. 
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Mass Diffusion Model of He Gas 

The He segregation process during freezing of liquid Li is delineated in Figure 7b. 
The governing equation, initial and boundary conditions for He gas mass diffusion 

process in the liquid region (x i S(t) ) are: 


9C +y dC = D ^C 


dt dx dx2 


(17) 


C(x, 0) = C 0 

(18a) 

J(6(t)) = (C(S,t)- Cl )V(t) 

(18b) 

^(X p t)= o 

dx 

(18c) 


Discretization of Equation (17) is done using Patankar’s high flux approximation 
[Patankar (1980)], namely: 

^cr 1 -cr> + j,.-j i ,..o. (19) 

The convection-diffusion flux at the left (Ji,-) and right ( Ji,+) boundary of node i are 
defined as: 

J i,sign = vCj + D^F(PeJ)(Cj - Cj +1 ) t (j = i for sign = + and j = i-l for sign 

where, 

D^ = 2D/(Ax j + Ax j+1 ) > 


F{PeJ) = <0,(l -0. l|P^| f> + < 0, - PeJ > f 

(20c) 

Pef = v/Df 
J J • 

(20d) 


= -) (20a) 

(20b) 
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The operator <a, b> is the largest of a and b. After linearizing the initial and boundary 
conditions in a consistent manner with that of the governing equation. He concentration 
in the liquid ahead of S-L interface, C(8,t), and the values of Q at each time step are 
determined by using a tri diagonal matrix solver. 

Volume Estimate of He Gas Voids 

As shown in Figure 7b, when the He gas concentration ahead of the S-L interface, 
C(5,t), reaches the HN concentration, He bubbles will nucleate and grow by absorbing 
gas atoms adjacent to them. The bubble growth will continue until all the gas atoms in 

excess of the solubility limit, or saturation concentration, Csat> have diffused into the 
bubbles. Because detailed analysis of the bubble nucleation and growth processes is 
beyond the scope of this work, it is assumed that He gas atoms in excess of Ch will be 
instantaneously contained in bubbles. The volume of these bubbles is calculated as: 



where. 


( 21 ) 


^ e = v( <C(X,t) * Ch)<k . (22) 

The numerical integration of the Equation (22) is carried out within the liquid region 
where the local He concentration, C(x,t), is greater than Ch. The program simulating the 
segregation process and the gas bubble evolution (SEGR.FOR) is listed in Appendix B-l 
and input data required to run the program (SEGR.INP) is listed in Appendix B-2. 

2.3.2 Results and Discussion 

By using the thermal properties of Li, the base case parameters listed in Table 3 
and the computer programs COPHASE.FOR and SEGR.FOR, die potential of He gas 
bubble formation by normal segregation and shrinkage of of Li during the freezing of Li 
coolant in the SP-100 primary coolant system is assessed and the results are presented 
and discussed in this subsection. 
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The segregation model described earlier is applied to three different coolant 
channel widths (0.01m, 0.05m, and 0.1m). The results are presented in Figures 8, 9 and 
10. Figure 8 delineates the calculated temperature profiles and Figure 9 shows the 
calculated He concentration profiles during the freezing process of liquid Li. It can be 
seen in Figure 8 that the temperature profile in the solid Li is almost linear. Also, the 
solidification velocities and the wall temperatures are independent of coolant channel 
width (see Figure 10). The solidification velocity and rate of change in the wall 
temperature are approximately 8 x 10-6 m/sec and -2 x 10-4 K/sec, respectively. Note 
that the freezing process of Li in SP-100 is extremely slow, suggesting a steady-state 
treatment of the problem could be sufficient Also, a constant wall heat flux 
approximation may be used in the analysis of freezing process of U with radiative heat 
transfer boundary. The low freezing velocity of Li is due to its high heat capacity and 
high latent heat of fusion as well as the low cooling rate at the wall. Figure 11 presents 
the calculated frozen fraction of the solidified liquid as a function of channel width, L. 

The results in Figures 12 and 13 show the volume of He gas bubbles formed by 
normal segregation and that of shrinkage voids during the freezing process of Li. The 
results indicated that the onset of He bubble nucleation occurs approximately 350 
seconds after the initiation of freezing at which the thickness of solidified Li is only about 
3 mm. These results suggest that formation of He gas bubble by normal segregation 
would occur throughout the PHTS of the SP-100 system during the Li freezing. It can be 
also seen in Figure 12 that the contribution of segregation bubble is only about 2% of the 
total volume of the void at the completion of solidification of each coolant channel. This 
corresponds to 0.055% of the volume of the coolant channels considered. 

One may argue that this amount of gas bubbles will not be sufficient to effectively 
remove the stress in the walls of PHTS during a subsequent restart of the SP-100. 
However, if these bubbles are captured by the freezing front into the solid Li, they will 
grow to compensate the volume shrinkage of Li. Thus, upon restart, the He 


22 



gas/shrinkage bubbles in the solid Li would effectively reduce the stress induced by the 
volume increase of Li. However, if the He gas bubbles are continuously driven by the 
solidification front, they will eventually collapse into the shrinkage voids. In order to 
minimize the stress during the restart of SP-100 and melting of Li, heating should be 
applied where the He/shrinkage void eventually form (see Figure 3). This would initially 
cause overheating of the wall (hot spot) due to the presence of the voids, but it will 
accommodate the increase in Li volume upon melting. It is not at all clear at present 
which of the two arrangement described in Figure 3 would occur in a solidifying Li in 
space. The subsequent research, therefore, focuses on the thermal effect of shrinkage 
void formation on the freeze and thaw processes of Li coolant in the SP-100 system. 
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(a) Temperature Profile 


Figure 8. Calculated Temperature Profile 
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Figure 10. Variation of Solidification Velocity and Wall Temperature 
with Time for Different Coolant Channel Widths 
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Figure 1 1 . Variation of Frozen Fraction for Different Coolant Channel Widths 
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Figure 12. 
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Evolution of Voids During Lithium Coolant Cooldown and Freezing 
(Coolant Channel Width = 0.01 m, Initial He Concentration = 

4.26 x 10 -6 kg/m3) 




Volume of Voids (m 3 ) 



Figure 13. Evolution of Voids During Lithium Coolant Cooldown and Freezing, 
(Coolant Channel Width = 0.05 m. Initial He Concentration 
= 4.26 x 10 -6 kg/m3) 
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2.4 Nomenclature 


English 

A w : Atomic weight (g/mole) 

BNP : Potential for bubble nucleation 

Cp : Heat capacity (J/kg.K) 

Ax : Spacial mesh size (m) 

C$* : Solubility limit in the solid (kg/m3) 

C Gas concentration concentration (kg/m3) 

D : Diffusion coefficient (m2/s) 

Ed : Activation energy for diffusion of a molecule of gas through the 

liquid (J) 

H : Enthalpy (J/kg) 

h : Boltzman constant (6.26 x 10 '34 j/s) 

hr : Radiative heat transfer coefficient (W/m2.K)W//m2.K) 

J : Nucleation rate in Equation (4) (Nucleations/s) 

J : Convection-diffusion flux (kg/m2.s) 

nt : Number density of gas molecules (molecules/m3) 

k : Thermal conductivity (J/m.K) 

L Coolant channel width (m) 

nne : Number of moles of He gas 

P : Pressure (Pa) 

Pe : Peclet number in Equation (20d) 

R : Gas Constant (8.3144 J/mole.K) 

icr : Critical radius of bubble (m) 

T : Temperature (K) 

t : Time (s) 
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V(t) 

V 

VHe 

v 

Xr 

x 

Greek 

6 

e 

♦ 

X 

e 

a 

r : 

p 

o 

o 

Superscripts 

e : 

n+1 

orm : 

nor 

ra -1 : 

Subscripts 
a 

SS : 

act 

e 


S-L interface velocity (m/s) 

Volume of the PHTS (m3) 

Volume of He bubbles (m3) 

Material velocity due to volume shrinkage (m/s) 
Location of shrinkage void/liquid (or solid) interface (m) 
Spacial coordinate (m) 

Solid crust thickness (m) 

Emissivity 

Dimensionless nodal temperature 
Latent heat of fusion (J/kg) 

Dimensionless nodal enthalpy or Wetting angle 
Frozen fraction of a node 

Ps/pL - 1 
Density (kg/m3) 

Surface tension in Equation (1) (N/m) 

Stefan-Boltzman constant (5.67 x 10 *8 W/ra^ K 4 ) 

Right boundary of a node 

Present time step 
Previous time step 

Ambient 

Supersaturation 

Actual 

Right boundary of a node 
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f 

HN, h: 


HTN 

ij 

o 

s 

L 

sat 

w 


At melting temperature 

Homogeneous nucleation 

Heterogeneous nucleation 

i-th and j-th node, respectively 

Initial condition 

Solid 

Liquid 

Saturation 

Left boundary of a node 
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Table 1. The SP-100 System Nominal Operating Conditions. 


Pressure 

Temperature 

Li Mass in the PHTS 

Li volume in the PHTS 

Gas Separator/ Accumulator volume 

Average He concentration 


3.5 x 104 p a 
1350 K 
110 kg 
255 Liter 
35 Liter 

4.26 x 10-6 kg/m3 


Table 2. Temperature Ranges and the Free Energies 
of Formation of NbT, NbN and LiT in SP-100. 


Reaction Temperature AF(kcal/mole) a 


(1) NbN Formation 

T>1273 K 

-50.1 

(2) Mixture of Nb2N 

673 <T< 1273 K 

-54.0 

and NbN Formation 



(3) LiT Formation 

773 <T <973 K 

- 16.7 

(4) Full Dissociation 



of LiT 

T > 973 K 


(5) NbT Formation 

T < 573 K 

-3.2 


a. The free energy of formation is the value at T = 298 K. 
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Table 3. Thermal Properties and Base Case Parameters 
Used in the Analysis of Helium Gas Segregation 
During Lithium Freezing 


Property 

Value 

liqyj&Lii 

Thermal conductivity 
Density 
Heat capacity 
Thermal diffusivity 

42.8 J/s.m.K 
517kg/m3 
4.384 x 103 J/kg.K 
1.89x10-5 m 2/ s 

■SplidLi 

Thermal conductivity 
Density 
Heat capacity 
Thermal diffusivity 

75.8 J/s.m.K 
530 kg/m3 
3.78 x 103 J/kg.K 
3.78 x 10-5 m 2/s 

Base Case 
Parameters 

Initial He concentration 
Diffusivity of He in liquid Li 
Diffusivity of He in solid Li 
Li Latent Heat of Fusion 
Li Freezing Temperature 

4.26 x 10-6 kg/m3 
5.96x10-9 m 2/s 
Negligible 
4.321 x 105 J/kg 
454 K 


34 




Table 4. Thermal Properties and Base Case Parameters used in the 
Analysis of Freeze and Thaw Processes of Li and LiF 


IHSHSSMI 

Liquid 

Solid 

Units 

Li k 

42.8 

75.8 

W/m.K 

P 

517 

530 

kg/m3 

Cp 

4384 

3780 

J/kg.K 

a 

1.89 x 10-5 

3.78 x 10-5 

m2/s 


Hsl 

4.321x105 

J/kg 

Tf 

454 


K 

LiF k 

1.30 

5.00 

W/m.K 

P 

1376 

1792 

kg/m3 

Cp 

2779 

2470 

J/kg.K 

a 

3.40 x 10-5 

2.94 x 10-5 

m2/s 

H S 1 

3.473 x 105 


J/kg 

Tf 

1122 


K 

PWC-11 k 


23.9 

W/m.K 

P 


6490 

kg/m3 

Cp 


323.4 

J/kg.K 

a 


1.14 x 10-5 

m2/s 

Tf 


2740 

K 

LiF vapor k 

0.017 


W/m.K 

He gas k 

0.1678 


W/m.K 

Base Case Parameters 

Ambient Temperature of Space 

250 

K 

Temperature of Wall or Heat Source 
Freeze: 

250 

K 

Thaw: 


1350 

K 

Thickness of PCM at Liquid State 

0.05 

m 

PWC-11 Wall Thickness 


7.6x10-4 

m 

Emissivity 

at outer surface of the wall 

0.8 


at inner surface of the wall 

0.5 



35 




3. THE EFFECTS OF SHRINKAGE VOID FORMATION ON THE FREEZE 
AND THAW PROCESSES OF LITHIUM COOLANT 


3.1 Physical Model 

As mentioned in the previous chapter, because the contribution of gas voids 
volume to total voids forms in the SP-100 system is small (-2%), its effect on the thermal 
behavior during the freeze and thaw processes is expected to be small and may be 
negligible. Therefore, in this chapter, only the effect of shrinkage void formation on the 
freeze and thaw processes is considered. Instead, cases with lithium-fluoride are 
analyzed for the comparison’s purpose and to satisfy the practical interests, since LiF 
exhibits much larger volume shrinkage upon freezing (23.2% compared to 2.5% for the 
case of Li) and it is known as an excellent phase change material (PCM) for thermal 
energy storage (TES) for solar dynamic power system. In addition, it is important to 
show the temperature variation in the heat transfer surface, the model being used in this 
analysis considers the effect of heat transfer wall. 

A simple numerical scheme is developed to simulate two different modes of 
shrinkage void formation in a finite one dimensional system, namely: (1) a shrinkage 
void forming at solid-wall interface (or wall void WV afterwards) (Figure 14a); and, (2) a 
shrinkage void forming at the center (or adiabatic boundary) of the system (center void or 
CV, afterwards) (Figure 14b). The former will occur in the absence of good adhesion 
between solidified PCM and the wall, while the latter is expected if liquid poorly wets the 
wall. Modeling this problem requires solving the heat transfer equations with two 
moving boundaries, namely, the solid-liquid (S-L) interface and the PCM-void (P-V) 
interface, subject to radiative heat transfer at the wall. 
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Because an analytical solution for such non-linear problems is not attainable, a 
numerical solution based on a single (solid/liquid) cell approach is used. The problem 
studied is that of a one dimensional slab initially filled with PCM at its fusion 
temperature in a container having a wall thickness, tc, symmetric about x = L (see Figure 
14). The wall is either radiatively cooled into space or radiatively heated from an 
isothermal heat source (e.g. a thaw heat pipe) at a temperature, T a . The other side of the 
slab, x = L, is thermally insulated. The heat transfer within the void, that is forming 
during the freezing process, is by conduction and radiation. The conductance of the void 
is assumed to vary inversely with the size of void, i.e. h V c = k v /8v [Wilson and Solomon, 
1986, Hunter and Kuttler, 1983]. The void is assumed to be filled with He gas and LiF 
vapor for the cases of Li and LiF, respectively. The sensible heat of the vapor is assumed 
negligible. Also, the temperature differential in the frozen crust is expected to be small 
because of the high thermal conductivity of Li. Hence, constant properties are assumed. 
During the freezing process, the location and velocity of the boundary “B” of the solid- 

represented by X r = yS(t), v(t) = y(dS(t)/dt), respectively, where y = (1 - Pl/p s ). In Figure 
14b, the boundary “B” representing the liquid- void interface moves in the opposite 
direction, hence Xr = L - T8(t), and v(t) = - r(d8(t)/dt), where r = (p$/pl - 1). 

During the thaw process, the location and velocity of the boundary “B” in Figure 
14a can be represented as, Xr = Xrf - y6(t), and v(t) = - y(d6(t)/dt), respectively, where, 
Xri is the initial thickness of void. In Figure 14b, X r = L - Xn + rS(t), and v(t) = 
r(d8(t)/dt). 

Governing Equations 

For the one-dimensional system shown in Figure 14 the temperature in the 
receding phases is uniform and equal to fusion temperature. The heat balance equations 
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in the developing phase (subscript d), and in the heat transfer wall (subscript c) are: 

Developing Phase 


Heat Transfer Wall 


PcCpc^ 


3T C \ 

C "dx/, 


Initial and Boundary Conditions 
T(x, 0) = T f , 

T s (S(t), t) = Tj_(S(t), t) = Tf , 


k s ^MS(t), t) - k L ^kS(t), t) = plhJ^ - v) 
dx dx ' dt 

and 


kc^-tct) =eo(T 4 (-t c ,t) - Tt) 
dx 


( 22 ) 

(23) 

(24) 

(25) 

(26) 
(27) 


forRgl4a _ (2g) 

^X/X_ *1 = 0 

3x , for Fig 14b. (29) 

In a recent work by Yang and El-Genk (1991), it has been shown that when a 
radiative heat transfer boundary condition is applied to the freezing of Li, a linear 
temperature profile develops within the solid crust and quasi-steady state approximations 
can be used. When the thickness of a developing phase is small, a linear temperature 
profile may be used. The governing equations are linearized using a fully implicit 
scheme as follows: 
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x=-tc 0 Xr S(t) L 

(a) Shrinkage Void Forming at Solid/Wall Interface 



x=-fc 0 S(t) Xr L 

(b)Shrinkage Void Forming at Center of System 


Figure 14. Schematics of Shrinkage Void Formation in a Solidifying Lithium 

Coolant in a SP-100 System 
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Developing Phase 


q'*qi 2 =[p d Cpc/l 2 + 1 6 n+1 - TjS 11 )- p d q 3r T f (8 n+1 - $ n )]/ At, 

where. 

(30) 

<{=-hoi(l? + 1 -T a ), 

(31) 


(32) 

%«*b 2 (Tf-T$ +1 ), 

(33) 

hot =(t? + TlJtlj 4 T a ) /[1/eo + 4T 1 /h 1 ], 

(34) 

hj-Sa 
tc ♦ 

(35) 

hl2 + 1/hy + IA 2 J f 

(36) 

*2 A 

0 * 

(37) 

h v -^.+ eo(l| + T4)(T3 + T4) 
Oy » 

(38) 

• T3 ' h 4£ + £) r, + *l 

(39) 


(40) 

Heat Transfer Wall 


<k -q^-P^pcO 1 At - Ii 

(41) 


The superscripts n and n+1 represent previous and present time steps. In order to 
solve for the temperatures of the wall and developing phase, the thickness of the 
developing phase, 5, should be known. An iterative procedure is used to determine 6 by 
solving the linearized heat balance at the S-L interface: 
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(42) 


^T, - T?* 1 ) « (Pa)p s H s i (g n ^- g°) > 

where, Pa is equal to 1 for freezing and -1 for thawing. 

The computer program developed to assess the effect of shrinkage void on the 
freezing process of Li and LiF (FREZ.FOR) and that for melting process (MELT.FOR) 
are listed in the Appendices C and D, respectively. 

32 Results and Discussion 

3.2. 1 Comparison with Available Analytical Solution 

The accuracy of the numerical scheme is verified by making a comparison with 
Wilson & Solomon’s analytical solution for a constant wall temperature (Wilson and 
Solomon, 1986). The thermal properties and base case parameters used in this 
comparison and in the subsequent analyses are shown in Table 4. In this comparison, it is 
assumed that conduction in the void is the only mode of heat transfer. The results are 
shown in Figure 15. This figure shows that the numerical solution is in reasonable 
agreement with the analytical solution for the base cases of interest (error of less than 
10%). However, the error in the freezing constant for the cases with a wall void is larger 
(as much as 26%) than that with a center void. Figure 15 also shows that for base case 
conditions (St= 1.78 and 6.20 for Li and LiF, respectively) freezing constants for the wall 
void cases are approximately 34% and 19% of those for the center void cases. It is also 
shown that in the center void cases, Li and LiF have similar freezing constants, while in 
the wall void cases, LiF has a much lower freezing constant than that for Li, due to the 
larger void size of the former. 

3.2.2 Freezing with Radiative Heat Transfer Boundary 

The results of the freezing process with a radiative boundary condition are shown 
in Figures (16) through (18). Figure 16 shows that because the difference between the 
fusion temperature of Li (454 K) and the ambient (250 K) is relatively small and the 
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thermal conductivity of U is high, the freezing velocity for Li is almost constant 
Conversely, because of the high fusion temperature of the LiF (1122 K) the temperature 
differential between the wall and ambient is approximately three times higher than in the 
case of Li, the freezing velocity of LiF is an order of magnitude higher. As the void size 
and crust thickness increase (both have lower thermal conductance than in the case of Li), 
the wall temperature decreases, resulting in a lower freezing velocity. This velocity is 
still an order of magnitude higher than that of Li (see Figure 16). This figure also shows 
that while the location of the void negligibly affects the freezing of the Li, it significantly 
influences the freezing velocity of LiF. As Figure 16 demonstrates, the freezing velocity 
for LiF, in case of a void forming at the wall, is about half that in the case where the void 
forms at the opposite adiabatic boundary. Figures 17 and 18 shows that because the LiF 
fusion temperature is much higher and void size is larger that those of Li, the ratio of 
conduction to radiation heat transfer through the LiF void is about four orders of 
magnitude lower than that for Li. 

3.2.3 Thaw with Radiative Heat Transfer Boundary 

The thaw process with a radiative boundary has a practical importance to space 
applications. The thaw process begins by assuming that a shrinkage void has formed 
either at wall or at the center of a system (or adiabatic boundary). The size of void is 
equal to that resulting from a complete solidification of liquid initially filled in a 
container, having a width L. In the base case, void thickness is 1.23 mm for Li and 11.6 
mm for LiF. 

The results are delineated in Figures 19 through 21. At the start of the thaw 
process of Li the heat transfer to the frozen crust is limited by the heat transfer through 
the void. This causes a large fraction of the heat input to be absorbed in the wall, 
resulting in a higher wall temperature (see Figure 19). However, as the Li melts, 
reducing the size of the void, a smaller fraction of the heat input is absorbed in the wall, 
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causing the rise rate of the wall temperature to decrease. Eventually, when the heat 
transferring through the void equals the heat input, the wall temperature reaches a 
maximum value (see Figure 20). Beyond this point, the wall temperature decreases with 
time as the void size becomes very small. Subsequently, the heat transfer rate through 
the void exceeds the heat input; the difference is compensated for by the decrease in the 
sensible heat of the wall. The variation of the heat flux ratio inside the void during the 
thaw process is delineated in Figure 21. 

In the case of LiF, because the initial void size is about an order of magnitude 
larger than that of Li, the thaw of LiF is always limited by the heat transfer through the 
void. This causes the wall temperature to rise monotonically as the thaw process 
continues. As shown in Figure 19, the maximum wall temperature during Li thaw (990 
K) is about two-thirds that of LiF (1336 K). The peak wall temperature during Li thaw is 
much lower than the melting temperature of the wall materials in the SP-100 (PWC-11), 
suggesting that a development of hot spots is unlikely during thaw in the SP-100 system. 
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Freezing Constant • 



Figure 15. Variation of Freezing Constant, X, with Stefan Number 
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Figure 17. Temperature Variation in Wall Void Cases During Freeze of Li and 
LiF (Radiative BC, Ta = 250 K) 
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Figure 19. Temperature Variation, with Radiative BC Wall Void, Ta = 1350 K, 
During Thaw Process of Li and LiF 
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Figure 20. Heat Flux Variation During Thaw Process of Li and LiF (Wall Void, 
Radiative BC, Ta = 1,350 K) 
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Figure 21. Variation of Heat Flux Ratio Inside Void During Thaw Process 
of Li and LiF (Wall Void, Radiative BC, Ta = 1350 K) 
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3.3 NOMENCLATURE 


English 

Cp 

H sl 

hr 


hvc 

k 

L 

st 

T 


Tl,T 2 : 


t 

tc 

V 

Xr 

Greek 

a 

r 

5 

8 » 

e 

P 

o 


Heat capacity (J/kg.K) 

Latent heat of fusion (J/kg) 

Radiative heat transfer coefficient (W/m2.K) 

Thermal conductance of the void (W/m2.K) 

Thermal conductivity (J/ra.K) 

Coolant channel width (m) 

Stefan number () 

Temperature (K) 

Temperatures at the center of the wall and developing phase, 

respectively 

Time (s) 

Thickness of the wall (m) 

Material velocity due to volume shrinkage (m/s) 

Location of shrinkage void/liquid (or solid) interface (m) 

Thermal diffusivity (m2/s) 

Ps/PL - 1 

Solid crust thickness, (m) 

Void thickness, (m) 

Emissivity 
Density (kg/m3) 

Stefan-Boltzman constant (5.67 x 10 -8 W/m2.K) 
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Subscripts 

c 

d 

f 

L 

s 

1 
v 
01 
12 
1 

2 

3 

4 


Wall region 

Developing phase 

Fusion 

Liquid 

Solid 

Initial 

Void 

Between wall and heat source (sink) 

Between the centers of wall and developing phase 

Within the half of the wall 

Within the half of the developing phase 

At the wall- void interface 

At the PCM-void interface 
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4. SUMMARY AND CONCLUSIONS 


It is shown that the formation of He gas voids either by HN or HTN during the 
cooldown process is unlikely. However, the formation of Nb-hydride and/or Nb-nitride 
in the SP-100 PHTS could increase the potential of Li-dewetting of the walls during 
cooldown, and hence the HTN of He bubbles. Analysis shows that HN of He bubble 
would occur by normal segregation during the freezing process of Li. Because of its high 
heat capacity and high latent heat of fusion as well as the low cooling rate of the wall, Li 
freezing is very slow. The results indicate that the total volume of He bubbles is insignif- 
icant (-2 %) compared to that due to volume shrinkage of Li during freezing. Hence the 
former is expected to play a minor role in relieving the stress in the walls during the re- 
start of the SP-100 system. Depending on where heat is applied relative to the location of 
the shrinkage void, stress may or may not be induced in the PHTS structure during the re- 
start of SP- 1 00 system. 

The effects of shrinkage void forming during freezing of lithium and lithium-fluo- 
ride on subsequent thaw processes are investigated using a numerical scheme that is 
based on a single (solid/liquid) cell approach. Results show that the formation or pres- 
ence of a shrinkage void at the wall-PCM interface substantially reduces the S-L interface 
velocity while voids, forming at the center of the coolant ducts, do not influence the 
freeze and thaw processes. 

Results also show that the presence of a void at die wall causes a substantial rise 
in the wall temperature during thaw. However, in the case of Li, the maximum wall tem- 
perature was much lower than the melting temperature of PWC-1 1, which is used as the 
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structure material in the SP-100 system. Hence, it is concluded that a formation of hot 
spots is unlikely during the startup or restart of the SP-100 system. 
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Appendix A-l 

List of Computer Program 
COPHASE.FOR 
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oooooooooooooooooooooooooooooooo 


1 


ft 


ft 

ft 


PROGRAM COPHASE 
IMPLICIT REAL'S (A-H. O-Z) 

COMMON /INPT1/ NT. TIME. NX. XL. HINIT , TINIT . TA. TP. TWALL. 

HSF. HLF. CONDSF. CONDLF. RHOSF . RHOLF. EMIT 
COMMON /CONVE1/ EPS1 . EPS2 . EPS3 

COMMON /OPT IN/ MLB, MRB. I LB. IRB. ICOMP. 1ENT, ITEM. SLOPT, 
PRTIME(20) 

COIAtON /MESHX/ XI (0:500). DX(500) . DT 

COLMON /DIMSV/ HX(0:500, 2). TX(0:500, 2). DELTA(0:2) 

COtAtON /OIMLV/ THETA(0:500, 2). PA!(0:500. 2). THEC(0:500) 
COMMON /LOCPR/ RHOX(0 : 500 , 2) , CPX(0 : 500 . 2) , CONDX(0: 500 . 2) 
COMMON /TEMPO/ STHETA(O:500). STHEC(0 : 500) . SPAI (0:500). 

SRHOX(0 : 500) , SCPX(0:500). SCONDX(0 : 500) . SHX(0:500). 
STX(0:500) 

COMMON /THERP/ RHO(50), CP(50), CONO(50). TEM(50). ENT(50) 
COOrfON/PROP/HL . HS .CL . CS . CONL . CONS .RHL.RHS.TFF 


NAMELIST /OPTION/ MLB. MRB. I LB, IRB. ICOMP. I ENT, ITEM, 
ft SLOPT. NPRT , NpStp 

*••• MLB - 1 FOR MOVING LEFT BOUNDARY 
MRB - 1 FOR MOVING RIGHT BOUNDARY 
I LB - INDEX FOR LEFT BOUNDARY CONDITIONS 

0 - ADIABATIC BOUNDARY 

1 - ISOTHERMAL WALL (TIME INDEPENDENT) 

2 - ISOFLUX WALL (TIME INDEPENDENT) 

5 - CONVECTIVE BOUNDARY 

4 - RADIATIVE BOUNDARY 
IRB - INDEX FOR LEFT BOUNDARY CONDITIONS 

0 - ADIABATIC BOUNDARY 

1 - ISOTHERMAL WALL (TIME INDEPENDENT) 

2 - ISOFLUX WALL (TIME INDEPENDENT) 

3 - CONVECTIVE BOUNDARY 

4 “ RADIATIVE BOUNDARY 

JCCCCC*** 9/26/1990 OPTION FOR IRB - 1 - 4 ARE NOT USED 


ICOMP - 1 FOR COMPARISON WITH ANALYTICAL RESULT IF AVAILABLE 
I ENT - 1 FOR UNIFORM INITIAL ENTHALPY DISTRIBUTION 

2 FOR NON-UNIFORM INITIAL ENTHALPY DISTRIBUTION 
(INPUT REQUIRED) 

ITEM - 1 FOR UNIFORM INITIAL TEMPERATURE DISTRIBUTION 

2 FOR NON-UNIFORM INITIAL TEMPERATURE DISTRIBUTION 
(INPUT REQUIRED) 

SLOPT - OPTION FOR LOCATION OF INTERFACE 

THREE DIFFERENT POSITION CAN BE CONSIDERED AS THE INTERFACE 

- 0. FOR DELTA BASED ON THETA(X.T) - 0 

- .5 FOR DELTA BASED ON THETA(X.T) - 1/2 

- 1. FOR DELTA BASED ON THETA(X.T) - 1 

PRTIME(I) - PRINT OUT RESULTS AT TIME PRTIME(I) (SEC) 

IF NO INPUT IS MADE FOR THIS VARIABLE. 

PRINT IS MADE ONCE AT THE TIME SPECIFIED BY 
INPUT VARIABLE 'TIME* IN THE NAMELIST * INPUT1 * . 

NAMELIST /INPUT1/ NT. TIME. NX. XL. HINIT. TINIT. TF. HSF. HLF. 

1 CONDSF , CONDLF. RHOSF. RHOLF. EMIT. TWALL. HWALL. TA. 

2 CPS CPL 
NAMELIST /CONVER/ EPS1 . EPS2, EPS3 
OPEN(10, ERR-9999. FILE-’ COPHASE. INP* , STATUS-’OLD* ) 

OPENfl 1 , ERR-9999. FILE-'COPHASE.VEL* , STATUS- * UNKNOWN • ) 
OPEN(12, ERR-9999. F I LE-* COPHASE. ENT* . STATUS- * UNKNOWN 1 ) 
0P£N(13, ERR-9999, FILE-'COPHASE.TEM* , STATUS- * UNKNOWN ’ ) 
0PEN(14. ERR-9999. FILE-’COPHASE.CRS* . STATUS- * UNKNOWN * ) 

C •••• DEFAULT VALUES FOR OPTIONS 

DATA ML8/0/, MRB/0/. ILB/1/. IRB/0/, SLOPT/O./. I FLAG/0/ 



C •••• CONSTANTS 

DATA PHI/3.1415926/, BOLTZ/5 . 670-08/ 

C •••• REAO INPUT 

READ(10. OPTION) 

READ! 10, INPUT 1 ) 

READ(10.CONVER) 

HL = HLF 
HS - HSF 
CL - CPL 
CS - CPS 
CONL - CONDLF 
CONS - CONDSF 
RHL - RHOLF 
RHS - RHOSF 
TFF - TF 

C •••• IDENTIFY TIME FOR PRINTOUT 
DPRT - T I ME/DF LOAT ( NPRT ) 

PRTIME(I) - 0. 

DO 5 1-1 .NPRT 

5 PRTIME( 1+1 ) - PRTIME(I)+OPRT 

C •••• MESH GENERATION: 

C NX - NUMBER OF NODES 

NX1 - NX - 1 
DX1 - XL/DFLOAT(NX) 

DX101 - DX1 • 0.1 

Xl(0) - 0. 

XI (NX) - XL 

DO 10 1-1, NX 
DX(I) - DX1 

10 XI(I) - Xl(I-l) + DX(1) 

C •••• THE TIME STEP SIZE: UNIFORM 
DT - TIME/DFLOAT(NT) 

C •••• PREPARATION FOR INITIAL CONDITIONS 

DELTA(O) - 0. 

DELTA! 1) - 0. 

DELTA(2) - 0. 

DDELTADT - 0. 

SDELTA - 0. 

HF - HLF - HSF 
HO - HWALL 

SGAM - 1 . - RHOLF/RHOSF 
IF (MLB . EQ. 1) SCAM - -SGAM 
IF (MLB . EO.0) SGAM - 0. 

C •••• INITIAL TEMPERATURE. HEAT CAPACITY. THERMAL CONDUCTIVITY. & 

C DENSITY 

C DIMENSIONLESS ENTHALPY. THETA(I.M). TEMPERATURE. PAI(I.M) 

C AND ENTHALPY CORRECTION FACTOR. THEC(I) 

DO 20 J - 1.2 
DO 20 I - 1. NX 
HX(I.J) - HINIT 

CALL PROPERTY(HX(I , J) ,TX(I , J) ,CPX(I , J) ,RHOX( I , J) ,CONDX(I , J)) 
THETA(I.J) - (HX(I.J) - HSFJ/HF 
THEC(I) - - THETA(I.J) 

IF(THEC(I) .GT. 0.) TH£C(I) - 0. 

IF(THEC(I) .LT. -1.) THEC(I) - -1 . 

PAl(l.J) - THETA(l.J) + THEC(l) 

20 CONTINUE 



o O o 60 6 606 C 4 Oo 


RH0IN=RH0X(1 ,1) 

WR1TE(*.1001) SNGL(RH0X(1 .1)) ,SNGL(RH0X(2 . 1 ) ) . 
t SNGL(RH0X(NX-1 , 1 )) , SNGL(RHOX(NX , 1 ) ) ... 

1001 FORMAT (/5X , ’DENSITIES’ , E12 . 5. 2X . E12 . 5 . 2X . E12 1 5i2X , E12 . 5) 


•••• INITIAL VALUES FOR CONSTANT TWALL< TF, PAIO - CONSTANT. 


IF(IL8 .EQ. 1) THEN 

CALL PROPERTY (HWALL , TWALL.CPXO .RHOXO .CONDXO) 

THETAO - (HWALL - HSF)/HF 

PAIO - CPXO» (TWALL - TF)/HF 

THECO - PAIO - THETAO 

HO - HWALL 

GOLD « 0. 

WRITE(*.21) SNGL(TWALL) .SNGL(THETAO) , SNCL(PAIO) 


*••• INITIAL VALUES FOR RADIATIVE BOUNDARY 


ELSEIF(ILB .EQ. 4) THEN 

THETOLD - THETA(I.I) 

TWOLD - TX(1 ,1) 

PAIOLD - PA1(1 .1) 

HROLD - EMIT«BOLTZ*(TINIT**2. + TA**2. )»(TINIT + TA) 
QO LD - HROLD*(TINIT - TA) 

STE - CPX( 1 , 1 ) *(TF - TA)/HF 
TA2 - TA»*2. 

WRITE(» ,22) SNGL(TWOLD) .SNGL(THETOLD) , SNGL(PAIOLO) 

END IF 


•••• TOTAL SYSTEM ENERGY [JOULE] 


ETOT - 0. 

DO 25 I - 1, NX 

ETOT - ETOT + RHOX(I .1 )*HX(I , 1 )*DX(I) 

ETOTOLD - ETOT 
ETOTINIT - ETOT 

QTOT - e. 

WRITE(12 .26) SNGL(TF) , SNGL(HF) .SNGL(HLF) .SNGL(HSF) . 

* SNGL(CPL) .SNGL(CPS) , 

4c SNGL(CONL) .SNGL(CONS) .SNGL(RHOLF) .SNGL(RHOSF) 

WRITE(» ,26) SNGL(TF) .SNGL(HF) .SNGL(HLF) .SNGL(HSF) , 

It SNGL(CPL) .SNGL(CPS) . 

* SNGL(CONL) .SNGL(CONS) .SNGL(RHOLF) .SNOL(RHOSF) 


•••• THE BEGINNING OF TIME LOOP ••••••••••••••••••••••••••••••••••< 

DO 9000 MA - 1 . NT 

M *>2 

C STORE PREVIOUS STEP VALUES IN THE TEMPORARY STORAGE VARIABLES 

DO 30 1-1. NX 
SHX(I) -HX(I.M-I) 

STHEC(I) - THEC(I) 

STHETA(I) - THETA(I.M-I) 

SPAI(I) - PAI(I.M-I) 

SRHOX(I) - RHOX(I.M-I) 

SCPX(I) -CPX(I.M-I) 

SCONDX(I) - CONOX(I.M-I) 

30 CONTINUE 

SDELTA - DELTA(M-I) 


"> o o 

• 

• 

• 

• 

OUTER ITERATION DELTA ITERATION 
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IDITER - 0 

DELTA(I) « DELTA(0) + DDELTADT*DT/2 . 
CONTINUE 

IDITER - IDITER +1 

DDELTA - DELTA(I) - DELTA(0) 

! ASSUMED CRUST THICKNESS 

c 



I 

! 

C 

c 

INNER ITERATION 


C 

C 

c 

c 

r 

CONVERGENCE ARE CHECKED IN EACH TIME STEP BY; 

SUMDIF - SUM(F(1.M)— SF(m«DX(I))**2.). I-ILEND. IREND 
SUMASS - SUM(SF(I)«DX(I))»*2. )), I-ILEND. IREND 

1 

I 

1 

1 

i 

V 

C 

c 

AND IF (SUMDI F/SUMASS)<- EPS1 . 
THE CONVERGENCE ARE ACHIEVED 


1 

1 

l 

V 

c 



I 
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NITER - 0 
SUMQINDT - QO LO 
CONTINUE 



NITER - NITER ♦ 1 




SUMDIF - 0. 
SUMASS - 0. 



c •••• 

r 

ADJUSTMENTS OF THE SIZE AND LOCATION OF 

THE BOUNDARY NODES 


c 

c •••• 

INITIAL LOCATION OF BOUNDARY NOOES 
IF (MLB .EO. 0 .AND. MRB .EQ. 0) THEN 
XLEND - 0. 

XREND - XL 
I LEND - 1 
IREND - NX 




c 

C •••• LOCATION AND SIZE OF MOVING LEFT BOUNDARY NODE; MLB - 1 
C 

ELSE IF (MLB .EQ. 1) THEN 

XLEND - SGAM*DELTA(M-1 ) 

XREND - XL 
I LEND - 1 

IF(XLEND .LT. XI(1)) THEN 
I LEND - 1 
GO TO 60 

ELSEIF(XLEND .GT. XI(NXI)) THEN 
I LEND - NX 
CO TO 60 


ELSE 

DO 50 IL-2.NX 

IF(XLEND .GT. XI(1L-1) .AND. XLEND -LE. XI(IL)) THEN 
I LEND - IL 
GO TO 60 

END IF 

50 CONTINUE 


oooo 


END IF 


60 DXI = XI (I LEND) - XLEND 

IF(DXI .LT. DX101 ) THEN 
I LEND *» IL +1 
DXI - XI (I LEND) - XLEND 

END IF 

DX(ILEND) « DXI 
C 

C •••• LOCATION AND SIZE OF MOVING RIGHT BOUNDARY NODE: MRB = 1 
C 

ELSEIF(MRB .EQ. 1) THEN 

XREND - XL - SGAM.DELTA(M-I) 

XLEND - 0. 

I LEND - 1 

IF(XREND .GT. XI(NXI)) THEN 
I REND - NX 
GO TO 62 

ELSEIF(XREND .LE. XI (1)) THEN 
I REND - 1 
GO TO 62 

ELSE 

DO 61 IL-1 ,NX1 

IF (XREND .GT. Xl(IL-l) .AND. XLEND .LE. XI(1L)) THEN 
IREND - IL 
GO TO 62 



ENDIF 

61 

CONTINUE 


ENDIF 

62 

0X1 - XREND - XI (IREND) 


IF(DXI .LT. DX101) THEN 


IREND - IL - 1 

DXI - XREND - XI (IREND) 

END IF 

DX( IREND) - DXI 

ENDIF 

K - I LEND 
L - IREND 

I LEND 1 - I LEND + 1 
1REND1 - IREND - 1 
IREND2 - IREND - 2 

DO 35 I - I LEND, IREND 

35 SUMASS - SUMASS + (HX( I ,M-1 )*DX( I) )»»2 . 


•••* PRODUCTION RUN: LEFT BOUNDARY NODE 

I LB - 1 FOR CONSTANT WALL TEMPERATURE AT LEFT BOUNDARY 

IF(ILB .EQ. 1) THEN 

C ALM - CPX(K— 1 ,M-1 )*DX(K— 1 )/CONDX(K— 1 ,M-1) 

ALO - CPX(K,M-1)*DX(K)/CONDX(K.M-1) 

ALP - CPX(K+1.M-1)*DX(K+1)/CONDX(K+1 .M-1) 

A1 - 2. *DT/(DX(K)*RHOX(K,M— 1 )) 

A2 - 1 ./(ALO + ALP) 


opoo oooo 


A3 - 1 ./ALO 
C A4 - 

c CHECK - STHETA(K) + A1 • (A2»PAI (K+1 .M-1 ) + A3.PAI0) 

c 

c I F (CHECK .GE. 0. .AND. CHECK .LE. 1.) THEN 

c THETA(K.M) - CHECK 

c THEC(K) - -THETA(K.M) 

c 

c ELSEIF(CHECK .GT. 1.) THEN 

c THEC(K) - -1 . 

c THETA(K.M) •» ( STHETA(K) + A1 • (A2»PA1 (K+1 .14-1 ) 

c ft + A3*PAIO - (A2+A3)»THEC(K)) ) / ( 1 .+A1 • (A2+A3) ) 

c 

c ELSE 

c THEC(K) - 0. 

c THETA(K.M) - ( STHETA(K) + A1 • (A2*PAI (K+1 ,M-1 ) 

c ft + A3»PAIO - (A2+A3)*THEC(K)) ) / (1 -+A1 »(A2+A3)) 

c ENDIF 

CHECK - STHETA(K)»srhox(k)/rhox(k ,m— 1 ) 

k -(1 .-8rhox(k)/rhox(k,m-1))«Hsf/Hf 

ft + A1*(A2*PAI (K+1 .14-1) + A3*PAIO) 

IF (CHECK .GE. 0. .AND. CHECK .LE. 1.) THEN 
THETA(K.M) - CHECK 
THEC(K) - — THETA(K.M) 

ELSEIF(CHECK .CT. 1.) THEN 
THEC(K) - -1 . 

THETA(K.M) - ( STHETA(K)»8rhox(k)/rhox(k.m-1 ) 
ft —(1 .— srhox(k)/rhox(k,m-1 ) )»Hsf/Hf 

ft + A1*(A2»PAI(K+1 ,M-1) 

ft + A3*PAIO - (A2+A3) *THEC(K) ) ) / ( 1 .+A1 . (A2+A3) ) 

ELSE 

THEC(K) - 0. 

THETA(K.M) - ( STHETA(K)*8rhox(k)/rhox(k.m-1 ) 
ft -(1 .— srhox(k)/rhox(k.m-1))*Haf/Hf 

ft + A1 • (A2*PAI (K+1 , M-1 ) 

ft + A3*PAIO - (A2+A3)*THEC(K)) ) / (1 .+A1.(A2+A3)) 

ENDIF 

HX(K.M) - HSF + THETA(K,M)*HF 

CALL PROPERTY(HX(K,M) ,TX(K,M) .CPX(K.M) , RHOX(K ,M) ,CONDX(K ,M) ) 
PAI(K.M) - THETA(K.M) + THEC(K) 

HD IF - ((HX(K.M) - HX(K.M-1))/HX(K.M-1))*»2. 


I LB - 2 FOR CONSTANT WALL HEAT FLUX BOUNDARY AT LEFT BOUNDARY 


ELSEIF(ILB .EQ. 2) THEN 
PRINT * 

PRINT CONSTANT WALL HEAT FLUX BOUNDARY * 


I LB - 4 FOR CONVECTIVETIVE BOUNDARY AT LEFT BOUNDARY 


ELSEIF(ILB .EQ. 3) THEN 
PRINT • 

PRINT CONVECTIVE BOUNDARY • 



oooo 


I LB - 4 FOR RADIATIVE BOUNDARY AT LEFT BOUNDARY 


ELSE I F( I LB .EQ. 4) THEN 

ALPA - BOLTZ»EMIT*DX(K)/CONDX(K ,M-1 )/2 . 

HRR - BOLTZ*EMIT*(TX(K,M-1 )**2. + TA2) • (TX(K .M-1 ) + TA) 
ft /(I. + 4.»ALPA*TX(K.M-1)*0. ) 

STE - CPX(K ,M-1 ) • (TF - TA)/HF 

C ALM - CPX(K— 1 ,M-1 )»DX(K— 1 )/CONDX(K— 1 ,14-1 ) 

ALO - CPXiK ,M-1 )»DX(K)/CONDX(K ,M-1 ) 

ALP - CPX(K4-1 ,M— 1 ) »DX(K4-1 )/CONDX(K4-1 ,M-1) 

A1 - 2. •DT/(DX(K)*RHOX(K,M-1 )) 

A2 - 1 ./(ALO + ALP) 

A3 - HRR/CPX(K,M-1)/2. 

C A4 = 

c CHECK - STHETA(K) + A1 *(A2«PAI (K+1 .U-1 ) - A3»St«) 

c 

c I F (CHECK .GE. 0. .AND. CHECK .LE. 1.) THEN 

c THETA(K.M) - CHECK 

c THEC(K) - -THETA(K.M) 

c 

c ELSEIF(CHECK .GT. 1.) THEN 

c THEC(K) - -1 . 

c THETA(K.M) - ( STHETA(K) + A1«(A2*PAI(K4-1 ,M-1 ) - (A2+A3) 

eft • THEC(K) - A3*Ste) )/(1 + A1*(A2+A3)) 

c ELSE 

c THEC(K) - 0. 

c THETA(K.M) - ( STHETA(K) + A1*(A2*PAI(K+1 ,M-1) - (A2+A3) 

eft • THEC(K) - A3*Ste) )/(1.+ A1*(A2+A3)) 

c END IF 

c 

CHECK - STHETA(K)»arhox(k)/rhox(k,m-1 ) 
ft —(1 .— 8rhox(k)/rhox(k,m-1 ))*Haf/Hf 

ft + A1*(A2*PAI(K+1 .M-1) - A3*Ste) 

I F(CHECK .GE. 0. .AND. CHECK .LE. 1.) THEN 
THETA(K.M) - CHECK 
THEC(K) - — THETA(K.M) 

ELSE IF (CHECK .GT. 1.) THEN 
THEC(K) - -1 . 

THETA(K.M) - ( STHETA(K)»«rhox(k)/rhox(k .m-1 ) 
ft -(1 .-arhox(k)/rhox(k,m-1 ))»Hsf/Hf 

ft 4- A1 *(A2*PAI (K+1 ,M-1 ) - (A2+A3) 

ft • THEC(K) - A3*Ste) )/(1.+ A1»(A2+A3)) 

ELSE 

THEC(K) - 0. 

THETA(K.M) - ( STHETA(K)*8rhox(k)/rhox(k,m-1 ) 
ft —(1 .-8rhox(k)/rhox(k.m-1 ))*Hsf/Hf 

ft 4- A1 *(A2»PA1 (K4-1 ,M-1 ) - (A24-A3) 

ft • THEC(K) - A3*Ste) )/(1.4- A1*(A24-A3)) 

ENDIF 

HX(K.M) - HSF 4- THETA(K,M)»HF 

CALL PROPERTY(HX(K,M) ,TX(K,M) .CPX(K.M) .RHOX(K.M) .CONDX(K.M) ) 
PAI (K, M) - THETA(K, M) 4- THEC(K) 


ENDIF I 

C OPTIONS FOR OTHER BOUNDARY CONDITIONS 



otiuo ouooooouo 


HD IF - ((HX(K,M) - HX(K,M-1))/HX(K,M-1))..2. 

SUMDIF « SUMD1F + ((HX(K.M) - HX(K .M-1 ) ) »DX(K) )*»2 . 
HRNEW - HR 


MAIN RUN: INTERIOR NODES 


DO 200 I *= I LEND 1 . I REND 1 

ALM - CPX( 1—1 ,M-1 ) »DX( 1—1 )/CONDX( 1—1 ,M-1 ) 

ALO - CPX(I ,M-1)*DX(I)/CONDX(I .M-1) 

ALP - CPX(I+1 ,M-1).DX(I+1)/CONDX(I+1 .M-1) 

A1 - 2.»DT/(DX(I)*RHOX(I ,M-1)) 

A2 - 1 ./(ALO + ALP) 

A3 - 1 ./(ALO + ALM) 

A4 - 

CHECK - STHETA(I) + A1 «(A2»PAI (1+1 ,M-1 ) + A3»PAI ( 1-1 .M) ) 

I F (CHECK .CE. 0. .AND. CHECK .LE. 1.) THEN 
THETA(I.M) - CHECK 
THEC(l) - — THETA(I ,M) 

ELSEIF (CHECK .GT. 1.) THEN 
THEC(I)«-1. 

THETA(I.M) - (STHETA(I) + A1 • (A2*PAI ( 1+1 .*A-1 ) + A3«PAI ( 1-1 .M) ) 
ft - A1.(A2+A3)*THEC(I))/(1. + A1«(A2+A3>) 

c ELSE 

c THEC(l) - 0. 

c THETA(l.M) - (STHETA(I) + A1 • (A2«PA1 (1+1 .M-1 ) + A3«PAI ( 1-1 .M)) 

c ft - A1«(A2+A3)*THEC(I))/(1. + A1*(A2+A3)) 

c ENDIF 

c 

CHECK - STHETA(I)*arhox(i)/rhox(i .m-1) 
ft —(1 srhox( I )/rhox( i ,ro-1))*Hsf/Hf 

ft + A1*(A2*PAI(I+1.M-1) + A3*PAI ( 1—1 ,M) ) 

IF(CHECK .GE. 0. .AND. CHECK .LE. 1.) THEN 
THETA(I.M) - CHECK 
THEC(I) --THETA(I.M) 

ELSEIF (CHECK .CT. 1.) THEN 
THEC(I) - -1. 

THETA(I.M) - (STHETA(I)*srhox( I )/ rhox( I .m-1 ) 
ft -(1 .— srhox( I )/rhox( i ,m-1))*Hsf/Hf 

ft + A1*(A2*PAI(I+1 .M-1) + A3*PAI(I-1.M)) 

ft - A1*(A2+A3)*THEC(I))/(1. + A1*(A2+A3)) 


ELSE 

THEC(I) - 0. 

THETA(l.M) - (STHETA(I)*srhox( I )/ rhox( i . ra — 1 ) 
ft -(1 .-arhox( i )/rhox( i ,m-1))»Hsf/Hf 

ft + A1*(A2*PAI(I+1 ,M-1 ) + A3*PAI(I— 1 ,M) ) 

ft - A1*(A2+A3)*THEC(I))/(1 . + A1*(A2+A3)) 

ENDIF 


HX(I.M) - HSF + THETA(I ,M)»HF 

CALL PROPERTY (HX( I ,M) ,TX(I ,M) ,CPX(I ,M) .RHOX(I ,M) ,CONDX(I ,M)) 
PAI(I.M) - THETA(I.M) + THEC(I) 



o o o 


C •*.. CHECK CONVERGENCE FOR THETA OR THEC WHICH EVER IS CONVENIENT 
HDIF - ((HX(I.M) - HX(l,M-1))/HX(I,M-1))..2. 

SUMDIF = SUMD1F + ((HX(I.M) - HX( 1 .W-1 ) ).DX( 1 ) ) . 

200 CONTINUE 

• ••• CLOSING RUN: RIGHT BOUNDARY NODE (IDIABATIC BC ONLY) 


DX01 - DX(L) + DX(L-I) 

ALM ** CPXfL-1 ,M-1).DX(L-1)/C0NDX(L-1 ,M-t) 

ALO « CPXfL.W-1)»DX(L)/CONDX(L.M-1) 

C ALP - CPX(L+1 ,M-1 )*DX(L+1 )/CONDX(L+1 ,M-1 ) 

A1 - 2.»DT/(DX(L)*RHOX(L,M-1 )) 

A2 - 1./(ALO + ALO) 

A3 - 1 ./(ALO + ALM) 

C A4 • 

c CHECK - STHETA(L) + A1 *A3*PAI (L-1 ,M) - A4*THETA(L— 1 ,M) 

c IF (CHECK .CE. 0. .AND. CHECK .LE. 1.) THEN 

c THETA (L.M) - CHECK 

c THEC(L) - -THETA(L.M) 

c ELSEIF(CHECK .CT. 1.) THEN 

c THEC(L) ■ -1 . 

c THETA(L.M) - (STHETA(L) + A1»A3»(PAI(L— 1 ,M) - THEC(L)) 

c A - A4*THETA(L-1.M))/(1. + A1«A3) 

c ELSE 

c THEC(L) - 0. 

c THETA(L.M) - (STHETA(L) + A1 «A3» (PAI (L-1 ,M) - THEC(L)) 

e * - A4*THETA(L-1,M))/(1. + A1*A3) 

c ENDIF 

CHECK • STHETA(L)*»rhox(L)/rhox(L,m-1) 
k -(I .-8rhox(L)/rhox(L,m-1))*Haf/Hf 

* + A1«A3»PAI(L— 1 ,M) - A4*THETA(L-1 ,M) 

IF(CHECK .CE. 0. .AND. CHECK .LE. 1.) THEN 
THETA(L.M) - CHECK 
THEC(L) - -THETA(L.M) 

ELSEIF(CHECK .CT. 1.) THEN 
THEC(L) ™ —1 . 

THETA(L.M) - (STHETA(L)«srhox(L)/rhox(L.m-1 ) 
k -(1 .-arhox(L)/rhox(L,m-1 ))»Hsf/Hf 

k + A1*A3»(PAI(L— 1 ,M) - THEC(L)) 

k - A4*THETA(L-1,M))/(1. + A1.A3) 

ELSE 

THEC(L) - 0. 

THETA(L.M) - (STHETA(L)«srhox(L)/rhox(L.m-1 ) 
k -(1 .-srhox(L)/rhox(L,m-1))»Haf/Hf 

k + A1 *A3* (PAI (L-t ,M) - THEC(L)) 

k ~ A4*THETA(L-1 ,M))/(1. + A1*A3) 

ENDIF 

HX(L.M) - HSF + THETA(L,M)*HF 

CALL PROPERTY(HX(L,M) ,TX(L,M) ,CPX(L,M) ,RHOX(L,M) .CONDX(L.M)) 


PAI (L, M) = THETA(L, M) + THEC(L) 


C *••• CHECK CONVERGENCE FOR THETA OR THEC WHICH EVER IS CONVENIENT 
HD IF = ((HX(L.M) - HX(L,M-1 ) )/HX(L,M-1 ) )**2 . 

SUMD1F - SUMD1F + ((HX(L.M) - HX(L.M-1 ) ) .DX(L) ) *»2 . 


C— 

C LOCATION OF INTERFACE 

C — — — 

C • ••• First, project THETA at x « 0. 
hi ■ dx(1) 
hi ** dxi2) 
h2 - dx(3) 
hihl - hi+hl 
hi h2 - hi + h2 
h12hi - hi + 2. •hi 
hi2h1h2 « hi + 2.*h1 + h2 
h22hi 2h1 - h2 + 2. •hi + 2. •hi 

theta(0,m) - theta(1 ,m)#(2. #h i ♦ hi )*h22h i 2h1/(h i hi *h i 2h1 h2) 
k - theta(2.m)#hi#h22hi2h1/(hih1*h1h2) 

k + theta(3.m)%hUh12hi/(hi2h1h2*h1h2) 

c If (delta(m-l) . le. 0. .and. theto(0,m) .ge. SLOPT) then 

c del ta(m) • 0. 

c go to 550 

c end if 

ISL1 - 1 

c i I endml ■ i I end - 1 

DO 541 I - I LEND. IREND1 

IF(THETA(I.M) .GE. 0. .AND. THETA(I.M) .LT. 1.) THEN 
ISL1 - I 

DELTA(M) - X1(1SL1) - DX( ISL1 ) •THETAflSLI ,M) 

C GO TO 551 

c ELSElf (THETA( I f M) .LE. SLOPT .AND. THETA( 1+1 ,M) .GT. SLOPT) THEN 

c IF(THETA(I,U) .LE. SLOPT .AND. THETA(I+1 ,M) .GT. SLOPT) THEN 

c ISL1 « I 

c DELTA(M) - XI(ISLI) - DX(ISL1)/2. + 2.*(SLOPT - THETA( ISL1 .M) )• 

c k (THETA( ISL1+1 ,M)-THETA(ISL1 .M) )/(DX( ISL1 )+DX( ISL1+1 ) ) 

c GO TO 551 

ENDIF 

541 CONTINUE 

551 CONTINUE 

I f (theta( l rend.rn) .le. elopt) then 
de I ta(m) ■» xrend 

print •, 'Complete solidification of the system* 
print • , *Ca I cu I at i on Terminated 1 
go to 9999 
c go to 545 

endi f 

c if(is!1 .eq . 1) then 

DELTA(M) - XI(ISLI) - DX(ISL1 )*THETA( ISL1 ,M) 

c ELSE 

c DELTA(M) - XinSLI-1) - DX(ISL1-1 )/2. 

c k + 2. •(SLOPT - THETA(ISLI-I.M))* 

c * (THETA(ISL1 ,M)-THETA(ISL1-1 .M) )/(DX(ISL1 )+DX(ISL1-1 ) ) 

c endlf 

C eumtheta ■ 0. 

C Do 555 I - I lend, istl 

C555 sumtheto ■ eumtheta + theto( i ,m)*dx( i ) 


c 

c 


ovgtheto « sumtheto/(x i ( i a 11 ) - xlend) 
delto(m) « (xi(iall) - x I end) *( 1 .-ovgtheto) 

c i a 1 1 « 1 

c DO 540 1 - 1 LEND , IRENDt 

c IF(THETA(I.M) .LE. SLOPT .AND. THETA( 1+1 ,M) .GT. SLOPT) THEN 

c 1SL1 » I 

c go to 545 

c END IF 

c540 CONTINUE 
c 

c545 xlmh - xi(iall) - dxfial1)/2. 
c xiph — xi(i*f0 + dx(iaM+1)/2. 

c 

c i f ( i a 1 1 .eq. 1) then 

c x*m3h » 0. 

c else 

c xim3h - xt(isll) - dx(iafl) - dx(ial1-1)/2. 

c endif 

c 

c if(ILB .eq. 1) then 

c upal - xrend 

c else 

c upal * 2. *x i ( i a 1 1 ) 

c endif 

c 

c unal - 6. 

c I i ter - 0 

c546 Xo - (upal+unal )/2. 
c liter- i lter+1 

c fx - theto(iaM-1 ,m)*(Xa - xirah)*(Xa - xiph) 

c k /((xim3h - ximh)*(xim3h - xiph)) 

c k +theto( I al 1 ,m)*(Xo - xim3h)*(Xa - xiph) 

c k /((ximh - x im3hW xlmh - xiph)) 

c k +theto( i a 1 1+1 ,mj* 1 X 0 - xim3h)*(Xa - ximh) 

c k /((xiph - xim3h)*(xiph - ximh)) 

c 

c wr i te(* ,2030) liter, lalt ,angl(Xo), angl(fx) 

c2030 format (2x, 'Niter-* , I5.3x, 'Nnode-' , I5,3x, 'Xa-\e12.5, 
c k 3x, 'fx-' t e12.5) 

c dfx - fx - SLOPT 

c if(doba(dfx) . le. eps3) then 

c delto(m) - xo 

c go to 550 

c elaelf(dfx .It. 0.) then 

c unai — xo 

c go to 546 

c e I ae 

c upal — xo 

c go to 546 

c endif 

c550 CONTINUE 

C WRITE(* ,609) SNG L ( THETA (K ,M) ) ,SNGL( THETA (K+1 ,M)) , 

C k SNGL(THETA(L-1 ,M)) ,SNGL(THETA(L,M)) 

CONVERGE - SUMDIF/SUMASS 

C 

C •••• CHECK ENERGY BALANCE 

C 

c 

C DETERMINE THE ENERGY OUTPUT WITHIN DT 

C FIRST, BY HEAT FLUX 


I F ( ] LB .EQ. 1) THEN 

ODXPAl = 2 . »(PAI (K,M) - PA10)/DX(K) 

QNEW = CONDX(K,M).HF/CPX(K,M)»DDXPAI 

ELSEIF(ILB.EQ.4) THEN 
QNEW « HRR.(TX(K,M) - TA) 

ENOIF 


SUMQINDT - SUMQINDT + QNEW 
GOUT - QNEW.DT 

C •••• THEN, BY TOTAL SYSTEM ENERGY [JOULE] 

ETOTNEW - 0. 

DO 705 I - I LEND , I REND 

705 ETOTNEW - ETOTNEW + RHOX( I ,M)»HX( I ,M)*DX( I ) 

ECOR - 0. 

DO 706 I - I LEND, IREND 

706 ECOR - (RHOX(I.M) - SRHOX(l)).HSF.DX(l) 

EOUT - ETOTOLD - ETOTNEW + ECOR 
EBALANCE - (QOUT - EOUT) / EOUT *100. 

EBALRAT - ( (QOUT - EOUT) / EOUT)*»2. 

TSTPRAT - OOUT/EOUT 

DIF - 0. 

DO 560 I - I LEND, IREND 

DIF - DIF + ((HX(I.M) - HX( I ,M— 1 ) )/HX( I ,M-1 ) )»»2 . 

560 CONTINUE 

IF(DIF .GT. EPS1 ) THEN 

DO 801 JI - I LEND, IREND 
HX(JI,M-1) - HX(J I ,M) 

CALL PROPERTY(HX( J I ,M-1 ) ,TX( JI ,M-1 ) ,CPX(JI ,M-1 ) . 
ft RHOX( J I ,M-1 ) ,CONDX( J I ,M-1 ) ) 

THETA(jl.M-l) - (HX( j I . M — 1 ) - HSF)/HF 
THEC(jl) - -THETA(jl.M-l) 

IF(THECQI) .GT. 0.) THEC(jl) - 0. 

IF(THEC(jI) .LT. -1.) THEC(jl) - -1 . 

PAI ( j I . M-1 ) - THETA( j I , U-1) + THEC(jl) 

801 CONTINUE 

GO TO 1050 

ENDIF 

WRITE(*. 610) SNGL(EOUT) , SNGL(QOUT) , SNGL( EBALANCE) . 
ft SNGL(TOTEOUT) , SNGL(QTOT) , 

ft SNGL(TOTBAL) . I 8 1 1 , sng I (de I to(m) ) 

C WRITE(* ,602) MA. IDITER, NITER. SNGL(CONVERGE) . SNGL(SUMASS) 

C OUTER ITERATION DELTA ITERATION 

C CONVERGENCE TEST FOR DELTA 

C 

C IF(DELTA(M) .GT. 1 .D-5 .AND. DELTA(M-1 ) .GT. 1 .D-5) THEN 

C DIFDELTA - ((DELTA(M) - DELTA(M-1 ) )/DELTA(M-1 ) ) ..2 . 

C ELSE 

C DIFDELTA - OELTA(M) - DELTA(M-1 ) 

C ENDIF 

C 

C DELTA(M-I) - (DELTA(M-1)+DELTA(M))/2. 

C I f (del to(m-1 ) .le. edelto) del to(m-1 )»sdel to 

C 
C 


IF(DIFDELTA .GT. EPS2) THEN 



C DO 616 1 « I LEND. IREND 

C HX( I ,M-1 ) « HX(1 ,M) 

C CALL PROPERTY(HX(I.M-1),TX(I,M-1),CPX(I.M-1), 

C & RHOX( I ,M-1 ) , CONDX( I ,M-1 ) ) 

C 

C THETA(I.U-I) - (HX(I.M-l) - HSF)/HF 

C THEC(I) * — THETA( I ,M-1 ) 

C IFfTHECm .GT. 0.) THEC(I) - 0. 

C IF(THEC(1) -LT. -1.) THEC(I) - -1 . 

C PAI ( I , M-1) - THETA( I , M-1 ) + THEC(I) 

C81 6 CONTINUE 
C GO TO 1025 

C ENDIF 

C 

C »••• TRANSFORMATION INTO DIMENSIONED VARIABLES 

C 

DO 700 I « 1 . NX 

HX(I.M) - THETA(l.M)*HF + HSF 

CALL PROPERTY (HX( 1 ,M) ,TX(I ,M) ,CPX(1 ,M) ,RHOX( 1 ,M) , CONDX ( I ,M)) 

700 CONTINUE 

TOTEOUT - ETOTINIT - ETOTNEW 
OTOT - QTOT + OOUT 
TOTBAL - QTOT/TOTEOUT 

C •••• STORE PRESENT STEP VALUES IN THE PREVIOUS STEP VARIABLES TO ADVANCE 
C NEXT TIME STEP 

DO 710 1-1, NX 
THETA(I.M-I) - THETA(I.M) 

PAI(I.M-I) -PAI(I.M) 

TX(l.M-l) - TX(I.M) 

HX(I.M-I) - HX(I.M) 

RHOX(I.U-I) -RHOX(I.M) 

CPX(I.M-I) -CPX(l.M) 

CONDX(I ,M-1 ) - CONDX(I.M) 

710 CONTINUE 

DELTA(M-2) - SDELTA 
DELTA(M-1 ) - DELTA(M) 

DDELTADT - (DELTA(M-I) - DELTA(M-2) )/DT 

GOLD - QNEW 

ETOTOLD - ETOTNEW 

FROFRAC - DELTA(M)AL 

TIMEA - DT • DFLOAT(MA) 

T1MEB - DT • DFL0AT(MA+1) 

mal - mo/NpStp * NpStp 

if(mo .It. 20 .or. mal .«q. mo) then 

del taxi - delta(m)/XL 

WRITE(1 1 , 1133) SNGL(TIMEA). SNGL(DDELTADT) , SNGL(Del taXL) . 

* SNGL(Tx(ILEND,m)) 

WRITE(14, 1 130) SNGL(TIMEA), SNGL(OELTA(M) ) 
end! f 

DO 770 NK - 1 ,NPRT 

IF(TIMEA. LE.PRTIME(NK) .AND. TIMES .GT.PRTIME(NK)) THEN 
WRITE(12.712) SNGLiPRTIME(NK)) 

WRITE! 13,712) SNGL! PRT IME(NK) ) 

WRITE(13,713) SNGL ( FROFRAC ) 

713 FORMAT (2X, ’FROZEN FRACTION - *.E12.5) 

WRITE(* ,712) SNCL(PRTIME(NK)) 

ITIME - ITIME + 1 

WRITE( 12, 610) SNGL(EOUT) , SNGL(OOUT), SNGL (EBA LANCE) 



4 . SNGL(TOTEOUT). SNGL(QTOT), 

4 SNGL(TOTBAL). lall, ang I (de 1 1 o (m) ) 

WR1TE(*. 610) SNGL(EOUT) , SNCL(QOUT). SNGL(EBA LANCE) 
4 , SNGL(TOTEOUT) , SNGL(QTOT), 

4 SNGL(TOTBAL) , lall. ang I (de I to(m) ) 

WRITE( 1 2 , 1110) 

1110 rORMAT(/5X. ’DIMENSIONLESS ENTHALPY DISTRIBUTION'/ 

4 4X , ' •/ 

4 5X. 'POSITION ENTHALPY') 

DO 8010 IJ • ilend, i rend 

8010 WRITE(12, 1120) SNGL(XI ( I J) ) ,SNGL(THETA( 1 J ,M) ) 

WRITE( 13, 1111) 

1111 FORMAT (/5X, 'TEMPERATURE DISTRIBUTION'/ 

4 4X.' '/ 

4 5X. 'POSITION TEMPERATURE (K)') 

DO 8011 IJ «• ilend, i rend 

8011 WRITE( 13, 1120) SNGL(XI ( I J) ) ,SNGL(TX( I J ,M) ) 

1120 FORMAT (5X,E12.5,3X,E12.S) 

1130 FORMAT (5X,E12.5,3X,E12.5,3X,E12.5) 

END IF 

770 CONTINUE 

9000 CONTINUE 




4 

4 

4 

4 

712 

1133 


9999 


5X , * QOUT - \E12.5.' [JOULES]'/ 

5X . ' QOUT/EOUT - ' . E 1 2 . 5// 

5X , ’ NODE FREEZING®’ .15/ 

5X, 'CRUST THICKNESS ®’.E12.5//) 

FORMAT (/5X. 'COMPUTED RESULTS AT TIME -'.F12.2.' [SEC]') 
FORMAT(2x.E12.5.2x.E12.5.2x.El2.5.2x. E12.5) 


STOP 

END 

SUBROUTINE PROPERTY(H.T.CP.RHO.CON) 

IMPLICIT REAL*8 (A-H.O-Z) 

COLMON/PROP/HL . HS . CP L . CPS . CON L . CONS . RHOL . RHOS . T F 

IF(H .GT. HL) THEN 

T - TF ♦ (H - HL)/CPL 

CP - CPL 

RHO - RHOL 

CON - CONL 

ELSEIF(H .GE. HS) THEN 
T - TF 

CP - CPL - (HL - H)*(CPL - CPS)/(HL— HS) 

CON « CONL - (HL - H)*(CONL - CONS)/(HL-HS) 

RHO - RHOL - (HL - H)»(RHOL - RHOS)/(HL-HS) 

ELSE 

T - TF + (H - HS)/CPS 

CP - CPS 

RHO = RHOS 

CON - CONS 

END IF 

RETURN 

END 
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UNIVERSITY OF NEW MEXICO 



SOPTION 
MRB = 1 . 

I LB = 0 , 

IRB = 1 . 


NpStp *= 5, 
IDX - 0, 
SEND 
SINPUT1 
Fo - .5. 

Ft - .5. 

Xsbl « 0.05, 
NT «= 200. 

TIME ■ 2. , 

dt - .0005. 

Fx - 1 .05. 

NX ~ 325. 

XL *=.05. 


vopt = 0. . 

Av ® 7.80D-06, 

Bv =0. . 

o 

< 

n 

© 

Di v «* 1 . , 
NPRT - 10 
SEND 
SCONVER 
epsl - 2.D-2 
SEND 

CO - 4.26D— 06 . 
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SOPTION 
I LB = 4. 

$END 
$INPUT1 
NT = 13, 
H1NIT-1 .13206, 
T1NIT - 454. . 
HSF-0. 699906, 
CPS-3.963D3 . 
CONDSF-71 .28, 
RHOSF-530. . 
EMIT-0.8. 
TWALL-100. . 
SEND 
$CONVER 
EPS1 - 1.0-30, 
SEND 


NPRT - 10, 


TIME - 7800. . 

TF-454. . 
HLF-1.132D6, 
CPL-4.384D3. 
CONDLF-42 . 80 . 
RHOLF— 51 7 . . 

HWALL-.6D6 


EPS2 - 1 .0-12, 


NpStp - 1 


NX-50. XL - .05. 
TA-250. . 


EPS3 - 1 .0-14 



o o 


dtprt = t ime/df loot(nprt ) 

DO 5 1 = 1 , npr t 

PRTIME(I) = df loot ( i )*dtprt 
5 CONTINUE 

• ••• MESH GENERATION: STABILITY OF SOLUTION 
find maximum velocity, Vmax 
if(idx .eq. 0) then 
Vmax «* V( 1 ) 
do 10 i = 1 . nt 
if(Vmax .le. V(i)) then 
Vmax = V( i ) 
endi f 
10 continue 

dxl = 2 . *Oc/Vmax/d i v 

20 sum ■* 1 

do 25 i=2,nx 

25 sum « sum + f x**df I oat ( i-1 ) 

dxo - XL/sum 

X i (0) - 0. 

do 30 I - 1 . Nx 

dx(i) ■ dxo* f x**df I oat ( i -1 ) 

30 Xi(i) - Xi ( i— 1 ) + dx(i) 

i f (dxl . I t . dx( 1 ) ) then 
nx = nx+1 
go to 20 
end i f 

e I se 1 f ( 1 dx . eq . 1 ) then 
dxs m xsbl/Df I oat (Nx) 

DO 35 i « 1 .Nx 
dx( i ) « dxs 

35 Xi(i) » Xi(i-I) + dx(i) 

else 

Vmax - V(1) 
do 1011 i ~ l.nt 
if(Vmax .le. V( I ) ) then 
Vmax - V(i) 
end if 

1011 continue 

dxl m 2. *Dc/Vmax/d i v 

2011 sum - 1 

do 2511 i=2,nx 

2511 sum - sum + f x**df l oat ( i-1 ) 

dxo - XL/sun 

X i (0) - 0. 

do 3011 i - 1 , Nx 

dx(i) ■■ dxo«f x**df loot ( i— 1 ) 

3011 Xi(i) - Xi(i-I) + dx(i) 

I f (dxl . 1 1 . dx( 1 ) ) then 
nx m nx+1 
go to 2011 
end if 


endl f 

dxend « dx(l)/20. 
Nxl « Nx— 1 



OOU 5> u o OOU 


Nxl 0 e Nx-10 


C • ••• PREPARATION FOR INITIAL CONDITIONS 

SGAM - RHOS/RHOL - 1 . 

IF(MRB . EQ . 0 ) SGAM « 0. 

C •••• INITIAL CONCENTRATION 

00 40 J - 1 .2 
DO 40 I - 1 . NX 
Cx( 1 , J) - Co 
40 CONTINUE 


• ••• INITIAL TOTAL MASS OF GAS IN THE SYSTEM 

Citot - 0. 

DO 50 I - 1 , NX 

I Citot - Citot + Cx( I ,1 )*DX( I ) 

Ctotold - Citot 

Co*XL 

•••• Bubble Volume Colculotion 
ExcHe » Co*XL/AwHe 
VolHetot - ExcHe*Rg*To/Po 

CAdmaea - 0. 

Cf lux « 0. 

DELTA - 0. 
lend « 0 
tr - 0. 
i t I me * 0 

•••• THE BEGINNING OF TIME LOOP ••••••••••*••••«**«#•. 

9000 Cont i nue 

c DO 9000 MA « 1 , NT 

I t Ime - i t ime 4 1 

mo - 1 

tr - tr 4 dt 

V(ma) ■ Av 4 Bv*tr 4 Cv*tr**2. 

C • -(1 .4sGam)*V(ma) 

DELTA - DELTA 4 v(mo)*DT 
XI iq - Xlo-(1 .4sGom)*De I to 

I f ( x I i q . le. dXend) then 
lend ** 1 
endi f 

print •,* * 

print *,* S-L Front - * .Sng I (de I to) . • [M] • 
print Liquid Length *» ' .Sngl (XI iq) , * [M]' 

I f (MRB .eq. 1) then 

Xend - XLo - eGam«DELTA 
Vf - (XLo-Xend)/XLo 

print • ,* Void Forming At Right Boundary* 
print Void Location is *, engl(Xend),* [M]' 

print ♦ Vol (Void)/Vol (Sys) »* .engl(Vf) 


else 



op 9°? 2 oooo 


Xend *= XLo 
Vf = 0. 

print No Void Formotion’ 

print Void Locotion is sng((Xend),‘ [Mj ' 

print Vot(Void)/Vol(Sys) ='.sngl(Vf) 

endi f 


•••• Determine Equivalent Interfoce Diffusion Coefficients and 
Patonckor High Flux Function 


DO 60 i « 1 ,Nx 

D i ( i ) - 2. *Dc/(DX( i )+DX( i+1)) 
Pe - C/D i ( i ) 

Fpe(i) - FP(Pe) 

Cont inue 


•••• determine coefficients for left boundary node 


DFL - Di (l)»FPe(1) 

Dml - Dc/DX( 1 ) 

Pel - G/Dml 
FPel - FP(Pet) 

DFL1 - Dml »FPe1 
AwCl ) - 0. 

Ae(l) - — Fo*DFL 

IFflLB .EQ. 1) THEN I J- - -Cs«[d(del to)/dt] 

Am(1) - DX ( 1 )/d t + Fo» (DFL + C) 

S(1) *= Dx( 1 )/dt»Cx( 1 ,M-1 ) — Cs»v(ma) 

4 + (1 .—Fa)* (DFL*Cx(2.m-1) - (G + DFL)»Cx(1 .m-1)) 

ELSEIF(ILB .EQ. 2) THEN I J- - beto»(Ci - Cs) 

beto - 1./f1./v(mo) - DX(1)/2./0c) 

Am(1) - DX( 1 )/dt + Fo*(DFL + G - beto) 

S(1) •» Dx( 1 )/dt»Cx(1 ,M-1 ) - Cs*beto 
k + (l.-Fo)e (DFL»Cx(2.m-1) + (beto - 0FL - C)*Cx(1 ,m-1)) 

ELSEIF(ILB .EQ. 3) THEN I J- - 0. 

Am(1 ) - OX ( 1 )/dt + Fa*(DFL + G) 

S(1) «* Dx(1 )/dt »Cx(1 ,M-1 ) 

k + (l.-Fo)e (DFL«Cx(2.m-1) - (0FL + G)*Cx(1 .m-1 )) 

ELSEIF(ILB .EQ. 4) THEN I J- - (Co - Ca).[d(del ta)/dt] 

Am(1) - DX(1 )/dt + Fa* (DFL + G) 

S(1) - Dx(1)/dt»Cx(1 .M-1) + (Co - Cs)*V(ma) 
k + (1 .—Fa)* (0FL*Cx(2.m-1) - (DFL + G)*Cx(1 ,m-1 )) 

ELSE I F( I LB .EQ. 5) THEN I J- - (Cl - CS)*[d(del ta)/dt] 

Am(1) =* 0X(1)/dt + Fa* (DFL + G - v(mo)) 

S(1) - Dx(1 )/dt*Cx(1 ,M-1 ) - Cs*V(mo) 

4 + ( 1 .— Fo)* (DFL*Cx(2.m-1) - (DFL + G - v(ma))*Cx(1 .m-1)) 

ELSE I C(0. T) - Ch 

Am(1) - DX(1 )/dt + Fo*(DFL + C + DFL1 ) 

S(1) - Dx(1)/dt*Cx(1 ,M-1) + (G + DFL1)*Co 

4 + ( 1 .— Fo)* (DFL*Cx(2,m— 1) - (DFL + G + 0FL1)*Cx(1 .m-1)) 

ENDIF 


•••• DETERMINE COEFFICIENTS FOR INTERIOR NODES 


OOP 


c- 


DO 70 I - 2. Nxl 
DFi - Di ( i )»FPe( i ) 

DFi 1 - Di ( i— 1 )»FPe( i-1 ) 

Awfil - — Fo* (DF i 1 + G) 

Ae( i ) « -Fo»DFi 

Am( i ) m DX(i)/dt — (Ae(i) +Aw(i)) 

S(i) « Dx(i)/dt.Cx(i .M-1) - (l.-Fo)/Fo»(Aef i).Cx(i+1,m-l) 
It — (Ao( i )+Aw( i ) ) »Cx( i ,m— 1 ) + Aw( i ) »Cx( i-1 ,m-1 ) ) 

70 CONTINUE 


C 

C *••• DETERMINE COEFFICIENTS FOR RIGHT BOUNDARY NODE 

C — 

DFR - Di(Nx).FPe(Nx) 

DFR1 - Di (Nxl )»FPe(Nx1 ) 

Aw(Nx) - — Fo»(DFR1 + C) 

Ae(Nx) » 0. 


IFflRB .EQ. 1) THEN l J+ - Co*G 

Am(Nx) - DX(Nx)/dt + FoOFRI 
S(Nx) - Dx(Nx)/dt*Cx(Nx,M-1) - Co»G 
ft + (l.-Fo) * ( (G+DFR1)*Cx(Nx1 ,m-1) 

ft — DFR1 *Cx(Nx ,m-1 )) 

ELSEIF(IRB .EQ. 2) THEN I Ci+1 - Ci 

Arn(Nx) - DX(Nx)/dt -Aw(Nx) 

S(Nx) - Dx(Nx)/dt*Cx(Nx,M-1 ) 
ft + (1 .— Fo)/Fo«Aw(Nx)*(Cx(Nx1 ,m— 1) — Cx(Nx,m-1)) 


ELSEIF( IRB .EQ. 3) THEN I Ci+1 - Co 

Am(Nx) - DX(Nx)/dt + Fa*DFR - Aw(Nx) 

S(Nx) - Dx(Nx)/dt*Cx(Nx,M-1) + DFR»Co 
ft + (1 .-Fo)*((DFR1 + G)*Cx(Nx1 ,m-1) 

ft - (DFR + G + DFR1)*Cx(Nx,m-1)) 

ELSEIF(IRB .EQ. 4) THEN I Ci+1 - 0. 

Ani(Nx) - DXfNx)/dt + Fo.DFR - Aw(Nx) 

S(Nx) » Dx(Nx)/dt*Cx(Nx,M-1) 
ft + (1 .-Fo)»((DFR1 + C)«Cx(Nx1 ,m-1) 

ft - (DFR + G + DFRI).Cx(Nx.m-l)) 

ELSE I J+ - CoG + DiFl(CI-Co) 

Am(Nx) - DXfNx)/dt + Fo»(DFR + DFR1 + G) 

S(Nx) - Dx(Nx)/dt»Cx(Nx,M-1) - Co»(G - DFR) 
ft + (1.-Fo)*((DFR1 + G)»Cx(Nx1 ,m-1) 

ft - (DFR + DFR1 + G)»Cx(Nx,m-1)) 


ENDIF 


•••• GET SOLUTION FOR PRESENT TIME STEP; TRIDIAGONAL MATRIX SOLVER 


CALL TDMA(#Cx, Aw, An. Ae. S, Nx) 

C •••• STORE CALCULATED VALUES IN THE PRESENT VARIABLES 
DO 100 1-1, NX 
Cx(I,M) — sCx(I) 

c Cx(l.ffi) - Cx(i,m)/(1. - (1. + eGam)*v(m)*dt/xl ) 

dx(i) - dx(i)»(1. - (1. + sGam)*v(m)*dt/xl ) 


o o oooo 


x i ( i ) *» xi ( i-1 ) + dx( i ) 

ICO CONTINUE 

Cmox = Cx( 1 ,m) 

Imax ■■ 1 
Do 299 i - 1 .Nx 
if(Cmox .LE. Cx(i t m)) then 
Cmox « Cx( i ,ro) 

Imax « I 
end i f 

299 continue 

CxCh • Cx( Imax ,m)/Ch 

WRITE( 1 1 , 6030) SNGL(TR). CxCh 

i tck - i t ime/Npstp*Nps tp 

i f ( > t c k .eq. itime .or. itime .It. 20 .or. iend .eq.l) then 
WRITE( 14, 6030) SNGL(TR) . Delta 
end i f 

pr i nt • , * * 

print • , *T inie«* , SNGL(TR) , * Cx(1)-\ Sng I (Cx( 1 .m) ) 

Print •,*Max. C Appears ot i ••.imax,' Cx(Imax) 
k Sng I (Cx( Imax ,m)) 




.... 


CHECK MASS BALANCE 

DETERMINE MASS FLUX AT LEFT BOUNDARY 

I F ( I LB .EQ. 1) THEN I J- - -Ca» [d(de I to)/dt ] 

FluxL •* — Ca»V(mo) 

ELSEIF(ILB .EQ. 2) THEN 1 J- - beto*(Ci - Ca) 

FluxL m beto* (Fo»Cx( 1 ,»n)+( 1 .— Fo)»Cx( 1 ,m-1 ) — Cs) 

ELSE! F( I LB .EQ. 3) THEN 1 J- - 0. 


FluxL - 0. 

ELSEIF(ILB .EQ. 4) THEN I J- - (Co - Ca) • [d(<Je I to)/dt ] 

FluxL • (Co - Ca)»V(mo) 

ELSEIF(ILB .EQ. 5) THEN I J- - (Ci - Cs)»[d(del to)/dt] 

FluxL » (Fo»Cx(1 ,m) + (1 .— Fo)*Cx(1 ,n>-1 ) — Cs)»v(mo) 

ELSE I C(0. T) - Ch 

FluxL - (C + DFL1 )»Co - DFL1 »(Fo»Cx(1 ,m) + (1. - Fo)»Cx(1 .m-1 )) 

ENDIF 


.... DETERMINE MASS FLUX AT RIGHT BOUNDARY 

IF(IRB .EQ. 1) THEN I J+ - Co»G 

FluxR «» Co*G 

ELSEIF(IRB .EQ. 2) THEN I Ci+1 = Ci 

FluxR — G»(Fo»Cx(Nx ,m) + ( 1 .— Fo)»Cx(Nx .m-1 ) ) 

ELSEIF(IRB .EQ. 3) THEN I Ci+1 - Co 

FluxR • (G + DFR)* (Fo»Cx(Nx ,m) + (1 .-Fo)»Cx(Nx,m-1 )) 
& - DFR*Co 

ELSEIF( IRB .EQ. 4) THEN I Ci+1 - 0. 

FluxR ■ (C + DFR)*(Fo«Cx(Nx,m) + (1 .— Fo)*Cx(Nx,m-1 )) 


ELSE 


I J+ - CoG + DiFI(Ci-Co) 


6 o 


200 


C •••• 


2118 

2116 

2115 


205 

c 

c 


2051 


FluxR *= Co*(G - DFR) 

-f (G + OFR) • ( Fo*Cx(Nx ,m) + ( 1 .-Fa) ♦Cx(Nx ,m-1 ) ) 

END I F 


Obtoin total mass at the time step 

Citot « 0. 

DO 200 I «* 1 .Nx 

Citot * Citot + Cx( i t m)+dx( i ) 

cont i nue 

ovgC i » C i tot/X I i q 

Admoss ■ (FluxL - FluxR)*dt 
CAdmass ■ CAdmoss + Admass 
Errorl » Citot - Admoss - Ctotold 
BALI - Er ror 1/Ctoto I d *100. 

Bubble Volume Calculation 
nnuc «* 0 
sumdx ■ 0. 

Cexcess *» 0. 

i f (Cx( imox ,m) .GT. Ch) then 
if(istop . eq . 1 ) then 
wri te(12 , 2115) sngl(tr) 
do 2118 i - 1 , nx.2 
cxco » cx( i ,m)/co 
xixliq = x i ( i )/x I i q/100. 
write(12. 2116) xixliq. cxco 
format r5x.d20. 10.2x.d20 . 10) 
f ormot (5x , * t ime *» • t f10.3 i * sec') 

print •, ’optional stop at first occurrence of Cx(imax.m) 

stop 

end i f 

DO 205 i • 1 . Nx 

if(Cx(i,m) .GE. Ch) then 

Cexcess - Cexcess + (Cx(i.m) - Ch)*dx(i) 

end! f 

cont i nue 

if (Cexcess .GT. CAdmass) then 
Aval 1C «* Cexcess - CAdmass 
Avo i 1C - Cexcess 
CAdmass ■ 0. 
do 2051 i - 1, Nx 

if(Cx(i ,m) .GE. Ch) then 

Cx( i t m) « Ch 

Sumdx m Sumdx + dx(i) 

Inuc ■» 1 
Nnuc ■* Nnuc + 1 
Nucnode(Nnuc) *= i 
end i f 

cont i nue 
endi f 

i f (Aval 1C .GT. 0.) then 
THeMole • Avoi IC/AwHe 
HeMole * THeMol e/Trnuc 
AddVol - HeMole*Rg*To/Po 
Vo I He - Vo l He + AddVol 


> Ch* 


o o 


if(votHe .GT. VolHetot) then 

AdjCx = (VoIHe - Vo I He tot ) *Po*AwHe/ (Rg*To) 

Vo I He « VolHetot 
Redist « AdjCx/Sumdx 
DO 206 i = 1 , Nnuc 

206 Cx(Nucnode( i ) ,m) = Cx(Nucnode ( i ) .m) 4 Red i st »dx(Nucnode ( i ) ) 

e I sc 

AdjCx = 0. 
end i f 

end i f 

end i f 

VolSh - Delto*sGam 
Vo I tot - VoIHe 4 VolSh 

Writer, 6080) eng J (AddVo I ) . snglfVoIHe), engl(VolSh). engl(Voltot) 
WRITER, 6020) SNGLtC I tot ) . SNGL(Admaee ) . SNGL(Ctot o Id), 

& SNGL(Errorl). SNGL(BALl) 

STORE PRESENT STEP VALUES IN THE PREVIOUS STEP VARIABLES TO ADVANCE 
NEXT TIME STEP 
DO 300 1-1 .NX 
Cx(I.M-l) - Cx(I.M) 

300 CONTINUE 

Ctotold - Citot - Cexcess 4 AdjCx 
Avo i 1C — 0. 

0 

c **»« TRANSFORMATION INTO DIMENSIONLESS VARIABLES 


TB - Tr 4 DT 

CxCh - Cx( Imox ,m)/Ch 

itck « i t ime/NpStp*NpStp 

i f ( » t c k . eq. itime .or. I time . eq. 1 
k .or. Inuc .EQ. 1 .or. lend .eq.l) then 

Wr J te(13, 6090) engl(TR). 8 ngl(VolHe). sngl(VolSh). engl(Voltot) 
end I f 

print • , * * 

print •.’Time-*. SNGL(TR) , • Cx(1)«\ Cx(l ,m) 

DO 500 U • 1 .NPRT 

IF(TR .LE. PRTIME(IJ) .AND. TB .GT . PRTIME( I J ) ) THEN 
WRITER. 6050) SNGL(Tr) 

WRITE! 12.6050) SNGL(Tr) 

WRITER 1 2 , 6040) 

IPFIog » 1 

ENDIF |. 

500 CONTINUE I 

lf(IPFIog .eq. 1) then 
do 400 ik - 1 , nx ,2 
CxCh«Cx( i k ,m)/Ch 
CxCo - Cx(ik,m)/Co 
XiXLiq - Xi(ik)/XLiq 

Arg - (x i ( i k)— dx( i k)/2. )/(2. ♦DSQRT (Dc*t r)) 

Anal - (Ca*(1 . - DERF(Arg) ) 4 Co*DERF(Arg) )/Co 
400 wr i te(12.21 19) XlXLlq. CxCo. Anol 1 

2119 format (5x,d20. 10.2x.d20. 10.2x.d20. 10) 

IPFIog - 0 
endi f 
Inuc — 0 


i f ( i end . ne . 1 ) then 
go to 9000 
end i f 

c9000 CONTINUE 


C 

C FORMAT STMTS 

C 


5010 FORMAT (5X ,D12 . 5 , 3X, D12 . 5) 

6001 formot(1x,'L »* ,F10.4, * m*) 

6010 FORMAT(/5X. * INITIAL CONCENTRATION DISTRIBUTION'/ 

4 5x.El2.5.2X,El2.5.2X,£12.5.2X.El2.5) 

6020 FORMAT (//5X , 'MASS BALANCE AT THIS TIME STEP*/ 

4 4X * ' */ 

* 5X, ' Total Mass - \E12.5.‘ [Kg/M**3]‘/ 

4 5X,* Accumulation - *,£12.5/ 

4 5x,' Initial Mass ■ *,£12.5/ 

4 4X # * */ 

4 5x.* Difference • '.E12.5*' [Kg/M**3]’/ 

4 5X, * % ERROR - \E12.5) 

6030 FORMAT(3X.E12.5,3X,D20,10) 


6040 FORMAT (/5x ,* He Concentration Distribution*/ 

4 4x , * */ 

4 5x . *Pos i t i on(M) Concentration [Kg/M* •3]'/ 

4 4x.« *) 


6050 FORMAT (/5X. ‘COMPUTED RESULTS AT TIME ■ \E12.5.' [SEC]*) 


6070 FORMAT ( 3X . E 1 2 . 5 , 3X . 020 . 1 0 ) 

6080 formot( 5x , ’ Inf ormat i on for Void'/ 

& 5x! ’Added He Vol. «‘.e12.5.‘ [M««3]’/ 

& 5x. ‘Totol He Vol . -‘.e12.5/ 

* 5x,'Vol. of Shrinkage Void »‘,E12.5/ 

& 5x, ‘Totol Void Volume -‘.E12.5//) 

6090 f©rmot(2x,e12.5.2x,e12.5,2x.e12.5,2x,e12.5) 


9999 STOP 
END 


c TRIDIAGONAL MATRIX SOLVER : oWi.X(i-l) + oPi.X(i) + oEi.X(i+1) • Si 


SUBROUTINE TDMA(X, oW.oP.aE.S. N) 
IMPLICIT DOUBLE PRECISION(A-H.O-Z) 
DIMENSION X(1), oW(1),oP(1),aE(1),S(1) 
DIMENSION P(1000).Q(10O0) 
c 

c ( aW(1) - 0. , oE(N) - 0. ) 



DO 10 i«2.N 

DENUM - oP( i )+oW( i )*P( 1-1 ) 
P(l) - -oE(i)/DENUM 



Q(') - ( S(i)-oW(i).Q(i-l) ) /DENUM 
CONTINUE 

X(N) ■= O(N) 

DO 20 i *=N — 1 .1,-1 

X(i) - P(i).X(i+1)+Q(i) 

CONTINUE 

RETURN 

END 
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Input Data for Computer 
Program SEGR.FOR 
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ft 

ft 

c 

c 

c 

c 

ft 


c 

c 

ft 

c •••• 


c •••• 

1 

2 


PROGRAM SEGRGATE 
IMPLICIT REAL»8 (A-H, O-Z) 

DIMENSION PRTIME(20) , Nucnode(1000) 

COMMON /CONVE1/ EPS1 . EPS2, EPS3 

COMMON /MESHX/ Xi (0:1000). dx(1000), sdx(1000). DT 
COLWON /CONCE/ Cx(0: 1000,2), Oi(0:1000), FPe(0:1000), 
sCx( 1000) 

CO»A*>N /VELOC/ v(0: 1000) , T(0:1000), vel(0:2), tm(0:2) 

COMMON /MATRIX/ Aw (1000), Am(1000), Ae(1000), S(1000) 
NAMELIST /1NPUT1/ NT. TIME. NX, XL, Co. Cs, Ch. Dc. vopt, 

Av, Bv, Cv, Div, nprt, fx, ft. Fo, dt, xsbl . Co 
NAMELIST /CONVER/ EPS1 , EPS2, EPS3 

NAMELIST /OPTION/ MLB. MRB. I LB, 1RB, Npstp, IDX, istop 

IRB — 0 ; Adiobotic Right Boundary 

• 1 ; Constont supply of gos ot initiol concentration 

with the speed of (sol i di f i cot i on+vol ume shrinkage) 

Function F(Pe): Potonckor High Flux 
FP(z) - DMAX1 (0.D0 , 1 .00 — 0.1D0*DABS(z))»»5.D0 
+ DMAX1 (0.D0, -z) 

OPEN(10, ERR-9999. FILE-’SEGR6. 1NP' . STATUS-'OLO' ) 

OPEN! 1 1 , ERR-9999. FlLE-'SEGR6.Cmt ' . STATUS- ’ UNKNOWN * ) 
OPEN(12, ERR-9999. FILE-’SEGR6.Cxi 1 . STATUS- ’ UNKNOWN • ) 

OPENf 13, ERR-9999, FILE-’SEGR6.Vol ’ . STATUS- * UNKNOWN • ) 
0PENM4, ERR-9999. FILE»'SEGR6.Del ’ . STATUS- 1 UNKNOWN ' ) 

OPEN! 15, ERR-9999, FILE-’SEGRe.Tsm’ , STATUS- ’ UNKNOWN 1 ) 

OPEN (20. ERR-9999. FILE-’vel .Tem\ STATUS- * o I d • ) 

DEFAULT VALUES FOR OPTIONS 

DATA MLB/0/, MRB/0/. I LB/1/. IRB/0/. Vopt/0/ 

DATA dt/ 5./ 

CONSTANTS 

DATA PHI/3. 1415926/, Dc/5.96D-09/. Rhos/530./. Rhol/517./ 
DATA Ch/1 .7370-03/, Co/4.260-06/. Cs/1 .590-08/, 

Csot/1 .590-08/ 

DATA Po/7330./, To/454. /. Rg/8.3143/. AwHe/4. 00260-03/ 

READ NAMELIST INPUT 
READ (10, OPT ION) 

READf 10, INPUT1 ) 

READ (10, CONVER) 

write(11 ,6001) XL 

wr i te( 12,6001 ) XL 

wr I tei 13,6001 ) XL 

wr i te(14,6001) XL 

wr i te( 15, 6001 ) XL 

THE TIME STEP SIZE: UNIFORM 

DTin - TIME/DFLOAT(NT) 

IF(VOPT .GT. 0.) THEN 
DO 1 1-1, NT 

READ(20.5010) T(I) ,V(I) 

ELSE 

DO 2 1-1 .NT 

T(I) - dFLOAT(I)»dt in 

V(l) - Av + Bv»t(i) + Cv»t(i)*«2. 

END IF 

XLo - XL 
Nxo « Nx 


c *••• 


IDENTIFY PRINTOUT TIME 
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c 


wr i te( 6 , 6002) 
write(7. 6003) 
c wr«tef 8 . 6004) 

write(9, 6005) 
wr 1 te( 10 , 6006) 
wr I tel 17 . 6007) 
wri teM6. 6008) 
wri te(19.6009) 
end i f 

I f ( i PCM .ne. 1) then I LiF calculation 
C LIF properties 
Tf- Tff 
Ti-Tif 
xKe~xKsf 
xKI-xKlf 
Cps*Cpsf 
rhos«rhosf 
rhol-rhol f 
Hsl«Hslf 
xKv»xKvf 
endif 

n-5 

nx »2 

i f ( i opt . eq . 1 ) 90 to 9000 


c if(lvoid .eq. 1 ) then 

c do » 0.00001 

c slo » do 

c else 

do - 0 . 

«lo « do 
c endif 

doe - do 

do 10 i- 1 ,n 
To(l)«Ti 
Te(i)-Ti 
10 continue 

HI » Hsl + Cps*Tf 
eps3o - eps3 
epso • epsl 
eps 12 » dsqrt (epsl ) 

Com - 1 . - rhol/rhos 
rGom « Gam/( 1 .-Gam) 

As m xKs/rhos/Cps 
Ac ■ xKc/rhoc/Cpc 
al ■ rhoc*Cpc*tc/dt 
o2 « rhos*Cps/dt 

hc 8 - 2. *xKc/tc 

Al - Emi tout»Bol tz/hcs 

SumLam* 0 . 

do 50 !t-1 # nt 

ind » I t/npr*npr - I t 
time m df loot ( i t ) •d t 
Erdo m 1 .030 
lf(lt .eq. 1 ) then 
dmox ■ do 4 2.*eps3 




20 


e I se 

dmax o do 4 2.*vei*dt 
end i f 
dm in « do 
iout «* 0 

cont i nue 
iout *» iout 4 1 

C Assume ds; new solid crust thickness 
ds - (dmox+dmi n)/2 . 
dvs » rGom*ds 
if( ivoid .eq. 1) then 
sis « ds 4 dvs 
e I se 

sis - ds 
end i f 

C Equivalent heat transfer coefficients 
C At the outer wall of the cladding , hoi 

Oen - 1 ./Emi tout/Bo l tz 4 4 . *Ts(1 )**3/hcs 
hoi - (Ts( 1 )**2+Ta**2)*(Ts(1 )+To)/Den 

C within the void 

hss - 2.*xKs/ds 
hcvs — xKv/dvs 

hr2s - Emltin*Boltz*(Ts(4)**2 4 Ts(3)* *2)* (Ts(4) + Ts(3)) 
hv — hcvs -f hr2s 
c hv m hcvs 

IF(iBC .eq. 1) THEN 

C Overall between T1 and T2 

C2 This considers heat transfer across the void 

if(ivoid .eq. 1) then 
h12 - 1 ./(I ‘/lies 4 1 ./hv + 1 ./hss) 
e I se 

C BC1 This does not consider heat transfer ocross the void 
h12 « 1./(1./hcs 4 1 ./hss) 
end i f 

C BC1 Radiative Boundary 
awfl) • 0. 

ap(1) » hoi 4 h12 4 at 

oe(1) * -hi 2 

s(1) • oNTo(1) 4 ho1*Ta 
ow(2} ^ -h12 

ap(2) ■» h12 4 hss 4 a2*ds 
oe(2) » 0. 

s(2) • rhos«Hl/dt«(ds-do) 4 o2*do*To(2) 4 hss*Tf 
call tdma(Tn t ow f cp,oe , s . nx) 

dT5 - Al *(Tn( 1 )**4.— To**4)/( 1 .44. *AI *Tn( 1 )»0) 

Tn(5) - Tn(1 )-dT5 

C BC1 This considers heot transfer across the void 
lf(ivoid .eq. 1) then 

Tn(4) - hi 2* (Tn(1 )/hss 4 (l./hcs 4 1 ./hv)*Tn(2) ) 

Tn(3) • h12*(Tn(2)/hcs 4 (l./hss 4 1 ./hv)*Tn( 1 ) ) 

else 

This does not consider heot transfer ocross the void 


C BC1 



Tn(4) = h12»(Tnf 1)/hss + 1 ./hcs»Tn(2)) 
Tn(3) » h12*(Tn(2)/hcs + 1 ,/hs 9 .Tn( 1 ) ) 
end i f 


ELSE 

C BC2 Constant wail temperature 


Tn( 

,i) 

- Ta 

Tnj 

5 

« Tn( 1 ) 

Tn< 


- Tn(1) 


if(ivoid .eq, 1) then 

C BC3 Constant wall temperature, with void 
h12 - 1 ./(I ./hss + 1 ./hv) 
e I se 

C BC2 Constant wall temperature, without void 
h12 — hss 
endi f 

Tn(2) ■ (a2*do*Tof2) + rhos*HI • (ds-do)/dt 
k + hss*Tf + h12*Tn( 1 > )/( hi 2 + hss + a2*ds) 

i f ( i voi d .eq. 1 ) then 

Tn(4) » (hv*Tn(1) + hss*Tn(2) )/(hss + hv) 

e I se 

Tn(4) » Tn(l) 
end i f 

END IF 

C Check energy botance ot the system 
C Heat dux at the left boundary of the system 

Qis » -ho1*(Tn(1)-Ta) 

C Heat flux at the left boundary of the PCM 

Qlp * -h12*(Tn(2)-Tn(1)) 

C For constant wall temperature only 
i f ( i BC .eq. 2) Qls - Qtp 

C Heat flux at the right boundary of the system, x**xLB 4 delta 
C which is equal to Heat flux at the right boundary of the solid PCM 

Qrs • -rhos*Hsl •(ds-do)/dt 

C Totoi heat accumulated in the system within dt 
c Fully implicit scheme 

EQp •* (Qlp - Qrs)*dt 
EQs (Qls - Qrs)*dt 

C Energy accumulated within the cladding 
{ f ( i BC .eq. 1) then 
Ec « tc*rhoc*Cpc* (Tn(l ) - To(1)) 
e I se 
Ec - 0. 
endi f 

C Energy accumulated within the solid crust 

Ep m rhos*Cps*(ds*Tn(2) - do*To(2) )-(ds-do)» rhos*H I 

C Total Energy accumulated in the system 
Ea - Ec + Ep 


Ers * (EQs — Ea)/Eo 



Erp » (EQp - Ep)/Ep 


REcEp « Ec/Ep 
vel « (ds-do)/dt 

if(it .eq. 1 .or. ind .eq. 0) then 
c wr i te(* ,6080) iout, ang I (Ec) , ang 1 (Ep) , sng I (REcEp) , ang l (ve I ) 

endif 

C Check the convergence 

Cx » 2.*xK3*(Tf-Tn(2))*dt/rhoa/Hsl 
i f (Cx .LT. 0.) Cx « 0. 

Cy • do**2 + 4.*Cx 
dec • (daqrt (Cy)+do)/2. 

Ero » (de - dsc)/ds 
AEro ■ daba(Ero) 

Erdn - doba(da - dec) 

i f ( I t .eq. 1 .or. ind .eq. 0) then 
c write(*,6020) engl(Erdo), angl(Erdn). eng I (Erp) # eng I (Ere) 

endi f 

if (Erdn . LE. epsl) go to 30 

if(dec .GT. da) then 

dm in « da 

Erdo - Erdn 

go to 20 

else 

dmox « da 
Erdo - Erdn 
go to 20 
endif 

30 continue 

if(]vold .eq. t) then 

Qr In - -£mitin*Boltz*(Tn(4)**4 - Tn(3)**4) 

Qcin » -hcva*(Tn(4)-Tn(3)) 

qrat - Qcin/Qrin 

endif 

Stef - Cpe*(Tf - Tn(5))/Hel 
FC - (ela-8lo)/dt*dsqrt (t ime/As) 

SumLom « SumLom + FC 

do ■ de 

elo - 8(8 

dos ■ da 

do 40 i«1 ,n 



40 continue 

if(iopt .ne. 1) then 

if(M .eq. 1 .or. ind .eq. 0) then 

print • / * 
print • / • 

print •/Sol. converged' 
oQle - dobe(Qla) 
oOlp » daba(Qlp) 
oQra * dabs(Grs) 

write(7, 6010) eng I ( t ime) ,engl (Tn(1 ) ) 
write(10, 6010) engl f t ime) ,angl (Tn(2) > 
c write(8, 6010) eng I ( t ime) , eng I lTn(3) ) 

write(9, 6010; engl (t ime) ,engl (Tn(4) ; 



ooonn 


c write(6, 6010) sng I ( t ime) . sng I f Tn(5) ) 

writefll, 6010) sngl ( t ime ) , sngl (FC) 
writef12, 6010) sng 1 1 1 ime ) . sng I ( ve I ) 
write(13, 6010) sngl (time) , sng I f Ers ; 
write{14, 6010) sngl ft ime) , sng I (ds) 
writeflS, 6010) sng I f t ime ) , sng I f q rot ) 
write(16, 6010) sng I f t ime ) , sng l f REcEp) 
write(17, 6010; sng I f t ime ) , sng I (oQI s) 
write(18, 6010) sng I f t ime ) , sng I (oQI p) 
write(19, 6010) sng I ( t ime; f sng I (oQrs) 
end i f 
endi f 

50 continue 

ovFC » SumLom/DFI oat (nt ) 

wri te(10, 8010) 

wr i tef 1 1 ,8011) Sng I (Stef) 

wri tef 12.8012 ) 

wr i tef 13,801 3) 

write(14,8014) 

wr i tef 15,8015) 

wr i tef 16.801 6) 

wr i tef 19 ,8017) 

go to 9999 


9000 continue 

This proc. calculates Average Freezing . v-btont, avgFC, for given time 
period 

in this proc. Wall Temp., Tw, is varied from 250 K 
to 450 K for Lithium ; TffLi) •* 454 K 

to 1120 K for Lithium Fluoride ; Tf(LiF) ** 1122 K 
if(iPCM .eq. 1) then 
Tu m 450. 
e I se 

Tu - 1120. 
end i f 
T I = 250. 

dTe - (Tu-TO/25. 

T I ow « T I - dTe 


c 

c 

c 

c 


c 


110 


do 160 i Tw=1 .25 

Tw * Tlow + dTe*Df I oat ( iTw) 

To - Tw 

print *,'Tw« '.sngl(Tw) 

if(ivoid .eq. 1) then 

do ■* 0.00001 

slo » do 

e I se 

do « 0. 

slo « do 

end i f 

dos » do 


do 110 i»1 . n 


To 

Ts 


(it 


*T i 
<Ti 


cont i nue 


HI • Hsl + Cps*T f 
eps3o « eps3 
epso - epsl 
eps12 « dsqrt(epsl) 
Gam • 1. — rhol/rhos 



rGom « Gam/(1.-Gom) 


As » xKs/rhos/Cps 
Ac - xKc/rhoc/Cpc 
ol ■* rhoc*Cpc* tc/dt 
a2 * rhos*Cps/dt 

hcs ■» 2 . *xKc/tc 

Al * Em i tout *Bo l t 2 /hcs 

SumLom=0. 

do 150 it-1. nt 

ind « i t/npr*npr — it 
t ime * df I oa t ( i t ) *dt 
Erdo - 1 .030 
1 f ( i t .eq . 1 ) then 
dmax — do 4 2.*eps3 
e I so 

dmax — do 4 2.*vel*dt 

endif 

dm i n — do 

iout ■» 0 

120 continue 

iout - iout 4 1 

C Assume ds; new solid crust thickness 
ds - (dmox+dml n)/2 . 
dvs * rGom*ds 
if( ivoid .eq. 1) then 
s I s « ds 4 dvs 
e I se 

sis - ds 
end i f 

C Equivalent heat transfer coefficients 
C At the outer wall of the cladding, hoi 

Oen « 1 ./Emi tout/Bol 1 2 4 4. *18(1 )**3/hcs 
hoi - (Ts(1 )**2+To*»2) • (Ts( 1 )+Ta)/Den 

0 within the vo i d 

hss « 2.*xKs/ds 
hcvs • xKv/dvs 

hr2s - Em 1 1 i n*Bol tz«(Ts(4)**2 4 Ts(3)^2) • (Ts(4) 4 Ts(3)) 
c hv » hcvs 4 hr2s 

hv — hcvs 

IF(IBC .eq. 1) THEN 
C Overall between T1 and T2 

C2 This considers heat transfer across the void 

if(ivoid .eq. 1) then 
h12 « 1./(l./hcs 4 1 ./hv 4 l./hss) 
e I se 

C BC1 This does not consider heat transfer across the void 
M2 • 1./(1./hcs 4 l./hss) 
endlf 

C BC1 Radiative Boundary 
owf 1 } m 0 . 

opt 1 1 « hoi 4 h12 4 al 
oe(1) - -M2 

s(1) » a1*To(l) 4 ho1*To 



aw( 2) «= — h12 

op(2) “ h12 + hss 4 a2*ds 
oe(2) « 0. 

s(2) «* rhos*HI/dt*(ds-do) 4 o2*do*To(2) 4 hss*Tf 
coll tdmo(Tn, ow, op, oe , s , nx) 

dT5 - AI*(Tn(1)**4.-To**4)/(l ,44.*AI*Tn(1)**3) 

Tn(5) - Tn( 1 )— dT5 

C BC1 This considers heot transfer across the void 
iffivoid .eq. 1) then 

Tni4) • h12»(Tn(l )/hss 4 (1 ./hcs 4 1 ./hv)*Tn(2)) 
Tn(3) - M2*(Tn(2)/hcs 4 (1-/hss 4 1 ./hv) *Tn( 1 ) ) 
e I se 

C BC1 This does not consider heat tronsfer across the void 
Tnf4l • hf 2*(Tnf l)/hss 4 1 ./bcs*Tn(2 ) ) 

Tn(3) ■ h12*(Tn(2)/hcs 4 1 ./hss*Tn( 1 ) ) 
end if 


ELSE 

C BC2 Constant woH temperature 


Tn< 

,i) 

- To 

Tnl 


- Tnf 

Tn| 

3 

- Tn( 


if(ivoid .eq. 1) then 

C BC3 Constant wall temperature, with void 
h12 - 1 ./(I ./hss 4 l./hv) 
else 

C BC2 Constant wall temperature, without void 
ht 2 « hss 
end i f 

Tn(2) * (o2*do*To(2) 4 rhos*HI *(ds-do)/dt 
* 4 hss»Tf 4 h12*Tn(1))/(h12 4 hss 4 o2*ds) 

iffivold .eq. 1) then 

Tn(4) m (hv*Tn(1) 4 hss*Tn(2))/(hss 4 hv) 
e I se 

Tn(4) - Tn( 1 ) 
endi f 

END IF 

C Check energy bolance of the system 
C Heat flux at the left boundary of the system 

Qls - -ho1*(Tn(1)-To) 

C Heat flux at the left boundary of the PCM 

Qlp - —hi 2* (Tn(2)— Tn(1 ) ) 

C For constant wall temperoture only 
i f ( i BC .eq. 2) Qls - Qlp 

C Heot flux at the right boundary of the system. x«xLB 4 delta 
C which is equal to Heot flux ot the right boundary of the solid PCM 

Qrs ■ -rhos*Hsl*(ds-do)/dt 

C Totol heat accumulated In the system within dt 
c Fully implicit scheme 



EQp «# (Qlp — Qrs)*dt 
EQs = (Qls - Qrs)*dt 

C Energy occumuloted within the clodding 
if(iBC . eq. 1) then 
Ec « tc*rhoc*Cpc*(Tn(l ) - To(1)) 
e I se 
Ec o e. 
endif 

C Energy occumuloted within the solid crust 

Ep *= rhos*Cps* (ds*Tn(2) - do*To(2) )-(ds-do)* rhos*HI 

C Totol Energy occumuloted in the system 
Eo • Ec + Ep 

Ers - (EQs - Eo)/Eo 
Erp - (EQp - Ep)/Ep 

REcEp • Ec/Ep 
vet * (ds-do)/dt 

C Check the convergence 

Cx ** 2.«xKs«(Tf— Tn(2))*dt/rhos/Hs I 
i f (Cx .IT. 6 .) Cx - 0. 

Cy *> do**2 + 4 . *Cx 
dsc • (dsqrt(Cy)+do)/2. 

Ero o (ds - dscl/ds 
AEro « dobs(Ero) 

Erdn • dobs(ds - dsc) 

if(it .eq. 1 .or. i nd .eq. 0) then 
c wr i te( • ,6020) sngl(Erdo). sngl(Erdn). sng I (Erp) , sng I (Ers) 

end i f 

if(£rdn .L£. epsl) go to 130 

if (dsc .GT. ds) then 

dm i n « ds 

Erdo » Erdn 

go to 120 

e I se 

dmox « ds 
Erdo m Erdn 
go to 120 
endi f 

130 continue 

Stef ** Cps*(Tf - Tn(5))/Hsl 
EC = (s I s-sl o)/dt*dsqrt (t ime/As) 
c FC ■ st s/2./dsqrt (As* t ime) 

SumLom - SumLom + FC 
c print • # *it-\it, • Fc - *,FC 

do ■ ds 
slo * sis 
dos ■ ds 
do 140 i-1 , n 



140 continue 


i f ( i t .eq. 1 .or. ind .eq. 0) then 



c pr i nt • . * * 

c pr i nt • , * ' 

c print •,*Sol. converged* 

end i 1 

150 continue 

ovgFC * Sum Lcm/DF loot (nt ) 

wri te(11 .6010) eng I (Stef) , sngl (ovgFC) 

160 continue 

wr i te(1 1 ,8020) 

9999 stop 

6001 format(2x, ’Generate o plot . ‘/2x, ' input data.*) 

C6002 format (2x, * " * , *T1 * . “*•) 

6003 format (2x, “'* , '12’ . •"*) 

C6004 format(2x, ‘ , *T3* 

6005 format (2x , *"* , *T4‘ . 

6006 format (2x, ••** . ‘T5* . ••'*) 

6007 f ormat (2x, **** , *CI add » ng* , * ,% * ) 

6008 f ormot (2x , * '" , *PCM-Voi d* . * " * ) 

6009 format (2x, , *S-L Intf\“") 

6010 format(2x»E12.5.2x,e14.7) 

6020 f ormot (2x . * Erdo®' »E12.5,2x . ‘Erdn®* ,E12.5,2x, # Erp«*,el2.5,2x, 

& * Ere*' .612.5) 

6080 f ormot (///2x, * lout®* , 15/ 

& 2x . * Ec®* ,e12 .5,2x . * Ep®* ,e12.5,2x. *Rot io®* ,e12.5. 
k 2x, *vel®* ,e12.5) 

8010 format (2x. 'end of doto. */2x, ’ t i t le text is *,* M *, ’Temperature * . 

* K* , *"* , ‘ . */2x, *go. ‘ ) 

8011 format(2x. 'end of data. */2x, * t i t le text is * . * " 1 . * Freez i ng* . 

& • Constant. St - 0 .E12. 5. **",‘. */2x, * go. * ) 

8012 formot(2x» *end of data . */2x , * 1 1 1 le text is 0 . 0,,0 , , Vsl. M/sec*. 

k *"*,'. */2x, ’go. * ) 

8013 f ormot(2x , ‘ end of dato . # /2x, * t i t le text is *,*"', ' Error in * , 
k ‘Energy Balance* , 00,0 . 0 . */2x, ‘go. *) 

8014 format(2x. ‘end of doto. ‘/2x. ‘title text Is ‘ . 0 00 0 . ‘Del to, M',**** 

* ,*. V2x,‘go.‘) 

8015 f ormat (2x, ‘end of data. */2x , ‘ t I t le text is ‘ .““.'Heat Flux*. 
k ‘ Ratio inside Void, Qcond/Qrad * , * M ' , * . '/2x , * go . * ) 

8016 format (2x. ‘end of data. */2x, 0 t I t le text is • ,**'*, 'Energy Ratio*, 
k ‘ in the System, dEc/dEs* ,““,*. */2x, *go. * ) 

8017 format (2x, ‘end of dato. */2x, * t i t le text is ‘.““,‘Heot Flux at*. 

k * Various Position, W/M**2* . ‘ '**.*. */2x . *go. * ) 

8020 format(2x, ‘end of data. ‘/2x. ‘t i t le text is '.*• ,, , 'Freezing* . 
k * Constant vs Stefan Number* . * H * ,*. */2x, *go. * ) 

end 


c 



c TRIDIAGONAL MATRIX SOLVER : oWi.X(i-l) + oPi.X(i) + oEi.X(i+l) - Si 


SUBROUTINE TDMA(X , aW.oP.oE.S. N) 

implicit double precision (o-h, 0 - 2 ) 
DIMENSION X(1), 0W(l).aP(1).oE(1).S(1) 
DIMENSION P(1000) .0(1000) 
c 

c ( aW( 1 ) - 0. . oE(N) - 0. ) 

c 



c 

DO 10 i-2.N 

DENUM - aP(i)+aW(i)*P(i-l) 

P(i) - -oE(i)/DENUM 
Q(i) - ( S(i)-0W(i).Q(i-1) ) /DENUM 
10 CONTINUE 

X(N) - 0(N) 
c 

DO 20 i»N-1.1 . -1 

X(i) = P(i)*X(i+1)+Q(i) 

20 CONTINUE 
c 

RETURN 

END 



program freeze 

C This program uses dichotomic search for S-L front 
C 

C This program calculates temperatures T(i). i-1,5 
C i-1 , outer surfoce of the wall 
C -2, middle of the solid crust 
C -3, inner surfoce of clodding 
C *4, PCM-Vold interface 
C -5, outer surface of cladding 
0 

implicit double precision (o-h, o-z) 
character^ odum 

dimension Tst (2,2000) .Qst (3.2000) . Timst(2000), Vst(2000) 
dimension To(5),Ts(5). Tn(5) 

dimension aw(l000). op(l000), oe(1000), s(1000) 
data To/250./, dt/.l/. nt/2000/ f Emit in/. 8/. Emitout/.8/. 
k Bo 1 1 z/5 . 67D-08/ , npr/10/, xSL/1 ./. iTpr/500/ 

C UF properties 

k Tff/1122./, Tif/1122./. xKsf/4.9997/, xKlf/1.3/, 
k Cpsf/2469.74/, rhosf/1792 ./. rhol f/1376./. He I f/347271 ./. 
C Thermal conductivity of LiF vapor 
k xKvf/0. 016954/, 

C U properties 

k T f/454 ./, T i/454./. xKs/75.8/, xKI/42.8/, 

k Cps/3780 ./» rhos/530./. rhol/517./. Hsl/432100./, 

C Thermal conductivity of He gas, at 200 F 
k xKv/0. 16781/, 

C Properties of cladding material 

k t c/0. 00076/, Cpc/323 . 4/, r hoc/6490./. xKc/23.885/. 
k epsl/1 - 0 — 1 2/ # eps2/1 .75/, eps3/0.01/, tcm/2740./ 

C Initial size of void, thickness of the system 

data iopt/1/,iPCM/2/, » void/1/, iBC/2/, gl/25000./ 


c 

c 


c 

c 

c 


openf6, err-9999, 
openi 7, err-9999, 
openi 8, err-9999, 
©pent 9, err-9999. 
openi 10, err-9999 , 
openi 1 1 , err-9999, 
open(12, err-9999, 
openi 1 3 , e r r-9999 , 
openi 1 4 , e r r-9999 , 
openi 15, err-9999, 
openi 16, err-9999, 
openi 17, err-9999, 
openM 8 , e r r-9999 , 
openi 19, err-9999, 
openi 20, err-9999, 
open(21 , err-9999, 
openi 22 , e r r-9999 , 
openi 23, err-9999, 
open(24, err-9999, 

wr i te(l 1 ,6001 ) 

IF( i opt .ne . 1 ) the 



i 1 e- 

•T1 .dot' , 


i le- 

'T2.dat ' , 


i te- 

'T3.dat * , 


i le- 

'T4. dot ' , 


i 1 e- 

'T5.dat ' , 


i 1 e- 

*FC. dot * , 


i le- 

'Ve 1 .dot * 



'Ers.dat* 


i le- 

•SL.dat* . 


i le» 

•Qrot.dat 


i 1 e- 

•Erct.dat 


i le- 

'Qcld.dot 


i le- 

'Qsol .dot 


i le- 

•Qdel .dot 


i le- 

•Ttot .dot 


i le- 

*Qtot .dot 


i le- 

'Tpro.dat 


i le- 

'Tpro.dot 


i le- 

•Tpro.dat 


status-' unknown * 1 
status-* unknown ' ) 
statue-' unknown * \ 
status-' unknown* i 
status-' unknown * ) 
status-* unknown* 1 
status-' unknown ' ) 
status—’ unknown ' 1 
status-' unknown* I 
, status-' unknown ' ) 
, status-* unknown* i 
, status-' unknown* ) 
, status-’ unknown ' ) 
, status-' unknown ' ) 
, status-' unknown* ) 
, status-* unknown* I 
, status-' unknown * ) 
, status-'unknown* ) 
.status-' unknown*) 


wr I te(7, 6001) 
wr 1 tei 12,6001 ) 
wrf tei13,6001) 
wr I tei 14,6001 ) 
wr I tei 15,6001 ) 
wr i tel 16,6001 ) 
write(17.6001) 
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program me 1 1 

C This progrom uses dichotomic seorch for S— L front 
C 

C This program calculates temperatures T(i), i = 1 ,5 
C i— 1, middle of the wall 
C -2 , middle of the liquid 
C -3, inner surface of cladding 
C -4, PCM-Void interface 
C -5* outer surface of clodding 
C 

C during a melting processes 

implicit double precision (o-h, 0 — 2 ) 
chorocter*8 adum 

dimension Tst (2, 2000) .Gst (3.2000) . Timst(2000). Vst(2000) 
dimension To(5).Ts(5), Tn(5} 

dimension aw(1000), op(1000), oe(1000), s(1000) 
data Ta/1350 ./, dt/1 ./, nt/20000/. Emitin/.S/, Emitout/.8/, 
4 Bo It z/5. 670-08/, npr/10/, xSL/0.05/, iTpr/500/ 

C Li F propert ies 

4 Tf f/1 122./ , Tif/1122./, xKsf/4.9997/, xKtf/1.3/, 

4 Cpsf/2469 .74/, rhosf/1 792./. rho I f/1376 ./. Hs I f/347271 ./, 

C Thermal conductivity of LiF vapor 
4 xKvf/0 .016954/, 

C Li propert i es 

4 Tf/454./, Ti/454./, xKs/75.8/, xKl/42.8/, 

4 Cps/3780,/, rhos/530./, rhol/517./, Hs 1/4321 00 ./. 

C Thermal conductivity of He gas, at 200 F 
4 xKv/0. 16781/, 

C Properties of cladding material 

4 tc/0. 00076/, Cpc/323.4/, rboc/6490./, xKc/23.885/, 

4 epsl/1 .0-12/, eps2/1 .75/, eps3/0.0l/. tcm/2740./ 

C Initial size of void, thickness of the system 

data 1 op t /0/, 1 PCM/2/, lvold/1/. 1BC/1/, ql/30000./ 
data dv i/0 . 0001/, imachine/l/ 
if(imachine .eq.1) then 

openf6, err-9999, f i le-*TTV. dot * , status-*unknownM 
open(7, err-9999, f i le-'TQR. dat * , status** unknown * ) 
e I se 

open(6, err-9999, f i I e-*Tl .dot * , status** unknown* } 
openi7, err-9999, f i I e-*T2 . dat * , status*' unknown * ) 
open(8, err-9999, f i I e-' T3 . dot * , status-* unknown ' ) 
open(9. err-9999, f i I e-' T4 . dot * , status-’ unknown * ) 
open( 10 ,e r r-9999 , f i le«* T5 .dot * , status-* unknown’ ) 
openM 1 , err-9999, f i I e-* dv . dat * , status-* unknown * 1 
openM2 , err-9999 , f i le-’ Ve I . dat * , status-* unknown ’ 1 
opent 13, err— 9999, f 1 le-* Ers . dat * , stotus— * unknown ’ S 
openM4, err-9999, f 1 le-*SL.dat * , status-* unknown* 1 
openM 5 ,e r r-9999 , f i I e»*0rat . dat * , status- * unknown ’ 1 
openM 6 , err-9999, f i I e-* Erot .dat * . status-* unknown* ) 
openM7, err-9999, f i le-*Qc Id. dat * , status-* unknown ’ ) 
open M 8, er r-9999 , f i le-*Qsol .dot * , status-* unknown* S 
open( 19 ,er r— 9999 , f i le-’Ode I .dot * , status—' unknown ' ) 
endi f 


1 


format (2x, *Gene rate 0 plot . */2x, * Input data.') 

f ormat (2x , * " * , *T1 ’ , * *' * ) 

f ormat (2x ,****, 'T2* ) 

format (2x , *•*', *T3* , *•'*) 

f ormat (2x, •*** ,*T4* ,*”*) 

format (2x, * *" , *T5* , **‘ ' ; 

format (2x, **” , ’Clodding* . *"*) 

format (2x , * H * , *PCM-Void* , * '** ) 

format (2x ,**'*, *S-L Intf*.**") 


if(lmachine .ne. 1) then 



wr i te(6, 1 ) 
wri te(1 1.1) 
wr i te( 12, 1 ) 
wr i tel 13, 1 ) 
wr i te(14, 1 ) 
wri te 15,1 

wri te(16. 1 J 
wr i te(17 , 1 ) 

wr i te( 6 .2) 
wr i te(7,3) 
wri te(8,4) 
wri te(9,5) 
wri tel 10,6) 
wr i tel 17,7) 
wri tel 18,8) 
wr i te( 1 9 , 9) 
end i f 

n-5 

nx=2 

if(iPCM .ne. 1) then I LiF colculotion 
C LiF properties 
Tf« Tff 
Ti-Tif 
xKs-xKsf 
xKI-xKlf 
Cps«Cpaf 
rhos«rhosf 
rhol«rhol f 
Hs l**Hs 1 f 
xKv»xKvf 
end i f 


Hs m Cpa*Tf 
HI - He I + Ha 
CpI - Hl/Tf 
eps3o * eps3 
epeo — ©pel 
epe12 « dsqrt(epsl) 

Gom « 1 . — rhol/rhos 
rGam » (l.-Gam)/Gam 

do - 0.000 
doe « do 

0 Estimate the thickness of PCM and void 

C based on the assumption that the system was originally filled with liquid 
C at its fusion temperature 

dv i « xSL*Gam 
xPCM m xSL — dvi 
dvo « dvi 

do 10 i»1 , n 
Tof i)-Ti 
Ts(i)-TI 
10 continue 

All - xKI/rhol/Cp I 
As - xKs/rhos/Cps 
Ac - xKc/rhoc/Cpc 
ol • rhoc*Cpc* tc/dt 



a2 « rhol*Cpl/dt 


hcs ® 2 . *xKc/tc 

A1 « Emi tout*Bol tz/hcs 

do 3000 »t«1. nt 

i nd ** i t/npr#npr - i t 
t ime • df I oo t ( i t)*dt 
Erdo - 1.030 
if (it .eq. 1) then 
dmox • do 4 3.*eps3 
else 

dmox • do + 3.*ve»*dt 

endtf 

dm in » do 

lout ■ 0 

100 continue 

lout » iout 4 1 

C Assume ds; new solid crust thickness 
ds *» (dmox4dmi n)/2 . 

C Equivoient heot transfer coefficients 
C At the outer wall of the cladding, hoi 

Den ■* 1 ./Emi tout/Bol iz 4 4. *Ts(l )**3/hcs 
hoi • (Ts(l)**24Ta**2)*(Ts(l)4To)/t ) en 

C within the void 

dvs » dvo - (ds-do)*Gam 
if (dvo .It. 1 .d-06) then 
| vo id ® 0 
dvs m 1 . d— 06 
end i f 

his - 2.*xKI/ds 
hcvs » xKv/dvs 

hr2s • Emi t in*Bol tx*(Ts(4)**2 4 Ts(3)**2)*(Ts(4) 4 Ts(3)) 
hv m hcvs 4 hr2s 

C Overall between T1 ond T2 

C2 This considers heot transfer across the void 

IFfiBC .ne. 2) THEN 
if(ivoid .eq. l)then 
h12 - 1./(1./hcs 4 1 ./hv 4 1 ./his) 
else 

Cl This does not consider heot transfer across the void 

M2 - 1 ./(I ,/hcs 4 1 ./his) 
endi f 

C BC1 Radiative Boundary 
owfl) ■ 0. 
oe(l) - -M2 
ifNBC .eq. 1) then 
ap(l) • h ®1 4 h12 4 ol 
s(1) - ol *To(t ) 4 ho1*To 
else 

ap(1) - M2 4 ol 
s(l) - o1*To(1 ) 4 ql 
endi f 

0 wf2) - -hi 2 

opi2) * h12 4 hts 4 o2*ds 
oe(2) - 0. 


C BC1 



s(2) «= rhol *Hl/dt* (ds-do) + o2*do*To(2) + hls*Tf 

C BC1 

coll tdma (Tn . ow , op , oe , s , nx) 
i f ( iBC . ©q . 1 ) then 

dT5 «* Al *(Tn( 1 )**4 .-To**4)/(1 .+4, «AI *Tn( t )**3) 

Tn(5) - Tn(1)-dT5 
e I se 

Tn(5) ■ Tn(1) + ql/hcs 
endif 

C BC1 This considers heat transfer across the void 
if(ivoid .eq. t) then 

Tnf4) • h12*(Tn(1 )/hl s *f (l./hcs + 1 ./hv)*Tn(2)) 
Tn(3) «■ h12* (Tn(2)/hcs + (l./hls + 1 ./hv)*Tn(l)) 
e I se 

C BC1 This does not consider heot transfer across the void 
Tn(4) - h12*(TnO)/hls + 1 ./hcs*Tn(2) ) 

Tn(3) • h12«(Tn(2)/hcs + 1 ./h I s*Tn( 1 ) ) 
endif 

ELSE 

C BC2 Constant wol I temperature, without void 
if(ivoid .eq. 0) then 
h12 - h I s 

C BC3 Constant wall temperature, with void 
e I se 

h12 - 1 ./(I ./His + 1 ./hv) 
endi f 

Tnfl) ~ Tc 

Tn(2) ■ (o2«do*Tof2) + rho I *H 1 *(ds-do)/dt 
* + h I s*T f + hl2*Tn( 1 ) )/(h12 ♦ hls + o2*ds) 

C BC2 Constant wall temperature. 



C BC2 wi thout vo i d 

i f ( i voi d . eq . 0) then 
Tn (4) - Tn(t) 
e l se 

C BC2 with void 

Tn(4) m (hv*Tn( 1 ) + h I s*Tn(2))/(h I s + hv) 
end i f 

ENDIF 

C Check energy balance of the system 
C Heot flux at the left boundary of the system 

Qls • -ho1*(Tn(l)-Ta) 

C Heat flux at the left boundary of the PCM 
Qlp - -h12*(Tn(2)-Tn(1)) 

C For constant wall temperature heat flux ot the outer surface of 
C cladding should be equal to that for inner surface of clladding 

iffiBC .eq . 2) Qls «= Qlp 
I f 1 1 BC .eq. 3) Qls - ql 


C Heat flux at the right boundary of the system, x»xLB + delta 
C which Is equal to Heat flux ot the right boundary of the solid PCM 



Qrs = -h I s* (Tf-Tn(2) ) — rho I *Hs I • (ds-do)/dt 
C for melting process, heot conducted ot the RHS of liquid 
C only by temperoture grodient 
Qrsl - -hi s*(T f-Tn(2)) 

C Total heot accumulated in the system within dt 
c Fully implicit scheme 

EQp « (Qlp - Qrs)*dt 
EQs « (Ol s - Qrs) *dt 

C Energy accumulated within the cladding (Temperature variation in the 
C cladding is considered only when the rodiative ond constant heat flux 
C boundary conditions) 

if(iBC .ne. 2) then 

Ec * tc* rhoc*Cpc* (Tn(1 ) - To(1)) 

C For constant wall temp Ec * 0. 
else 
Ec - 0, 
endi f 

C Energy accumulated within the liquid PCM 

Ep - rho I *Cp I • (ds*Tn(2) - do*To(2) )-(ds-do)* rho I *Hs 

C Total Energy accumulated in the system 
Eo • Ec + Ep 

Ers - (EQs - Ea)/Ea 
Erp - (EQp - Ep)/Ep 

REcEp - Ec/Ep 
vel « (ds-do)/dt 

c if (it . eq. 1 .or. ind .eq. 0) then 

c wr t te(* ,6080) i t , iout ,sngl (Ec) .sngt (Ep) .sngl (REcEp) .sngt (vel ) 

c endif 

C6080 format (///2x , * 1 1=* , i5, * lout-* . 15/ 

c k 2x, *Ec-' ,e12.5,2x, 'Ep- # t e12. 5, 2x, *Rat io-* ,e12. 5. 
c k 2x , ' ve I-' ,e12. 5) 

C Check the convergence 

Cx « 2.*xKI*(Tn(2)-Tf)*dt/rhol/Hsl 
i f (Cx .LT. 0.) Cx = 0. 

Cy « do**2 + 4 . *Cx 
dsc - (dsqrt(Cy)+do)/2. 

Ero - (ds - dscl/ds 
AEro * dobs(Ero) 

Erdn * dabs(ds - dsc) 

c i f ( 1 1 .eq. 1 .or. ind .eq. 0) then 

c wr I te(* ,6020) sngl(Erdo), sngl(Erdn), sng I (Erp) , sng I (Ers) 

c endif 

c6020 f orroot (2x, # Erdo-\E12.5,2x, # Erdn«\E12.5,2x, # Erp-\e12.5.2x, 

c k 'Ers**' , e12 .5) 

If (Erdn .LE. epsl) go to 2000 

if(dsc .GT. ds) then 

dm in * ds 

Erdo • Erdn 

go to 100 

else 

dmax * ds 
Erdo *■ Erdn 
go to 100 



o u 


end i f 


2000 continue 

i f ( i vo i d -eq. 1) then 

Qrin *= —Em i t i n*Bo I t 2 * (Tn(4) **4 - Tn(3)**4) 
Qcin = -hcvs* (Tn(4)-Tn(3) ) 
qrat = Qcin/Qrin 
endi f 

do *» ds 
doe « ds 
dvo • dvs 
do 170 i«=1 , n 



170 continue 

if(it .le. 50 .or. i nd .eq. 0 .or. do .gt. xSL) then 

c print • ,* * 

print *, ' • 

print •. 'Sol . converged* 

oQle - dabs(Qls) 
oQlp - dobs(Qlp) 

c aQrs * dabs(Qrs) 

C for melting process, heot conducted ot the RHS of liquid 
C only by temperoture gradient 
oQrs *» dobs(Qrsl) 
if(imochine .eq. 1) then 

errite(6, 6200) sngl (time) ,«ngl (Tn(l)),sngl(Tn(2)) , sngl (Tn(4)) , 
k sngl (vel ) 

write(7. 6200) sngl ( t ime) , sngl (oQI s) , sngl (oQI p) , sng 1 (oQrs) , 
k sngl(Qrot) 

e I se 

write(7, 6010) sng I ( t Ime) , sng I (Tn(1 )) 
wri tef 10, 6010} sng I c t ime) , sng I f Inf 2) ) 
wrlteC8, 6010) sng I ( t ime) , sng I (Tn(3) ) 
write(9, 6010) sng 1 1 1 ime) , sng I (Tn( 4 n 
errite(6, 6010) sng I ( t ime j , sng l(Tn(5) ) 

writefll, 6010) sng I ft ime) , sng I (dvs) 
wrifce<12, 6010) sng let Ime) , sngl Cve ( 1 
err i te( 13 , 6010) sngl (time) .sngl(Ers) 
errite(14. 6010) sngl (time) .snglids) 
wrlte(15, 6010) sng I (t Ime) , sngl (qrat ) 
write(16, 6010) sng I ( t ime) , sng I (REcEp) 
write(17, 6010) sngl ( t (me ) ,sng I (aQI s) 
wrlteMB, 6010) sng I ( t Ime ) , sng If aQl p) 
wrfte(19, 6010) sng I ( t ime) , sng 1 (oQrs) 
endi f 
endi f 

6010 formot (2x , El 2.5, 2x, el 4. 7) 

6200 format (lx, f8.2,1x,4e14.$) 

lf(do .gt. xSL)then 
err i te(* ,2990) 

2990 format (//2x, 'Co I cu lot ion Terminated on Complete melting of, 
k 0 the System') 

go to 3001 
endi f 

if(Tn(1) .GE. tern) then 
wri te(e,299l) 

2991 format (//2x, 'Ca I cut ot ion Terminated on Wall Temp. Exceeds*. 
k 'Melting Temp. ') 

go to 3001 



end i f 

3000 cont i nue 

3001 continue 

if(imochine . ne . 1) then 

wr i te(10 , 8010) 

8010 f ormot (2x , * end of data . */2x . ’title text is * . * "' f * Tempe rot ure ’ . 

k \ K' , . */2x , * go. * ) 

wr i te(1 1 ,8017) 

8017 format (2x , ’end of data. */2x, 'tit le text is * , ' " * . 'Void* , 
k •; Size; M' '/2x ,' go. * ) 

wr i te(12,801 1 ) 

8011 format (2x. 'end of doto . '/2x , ’ t i t I e text is '."".'Vsl, M/sec* , 
k ' . '/2x. 'go. ' ) 

wr i te(13,8012) 

8012 f ormat(2x, 'end of data. */2x , ' t » t le text is '."".'Error in *. 
k 'Energy Bo I once '/2x ,' go ) 

wr i te( 1 4 , 8013) 

8013 f ormot (2x ,' end of data. '/2x . ' t i 1 1 e text is ’."",'Delto, M'."" 

* V2x.'go.') 

wri te(15,8014) 

8014 format(2x, 'end of dato . */2x, 'title text is '."".'Heat Flux', 
k ' Ratio inside Void. Qcond/Qrod '/2x, 'go. ' ) 

wri te(16,8015) 

8015 formot(2x, 'end of dato . '/2x , • t i t I e text is '."".'Energy Ratio'. 

It ' in the System, dEc/dEs '/2x. 'go. ' ) 

wr i te( 19 ,8016) 

8016 formot(2x, 'end of data. */2x, ’ t i t I e text is '."".'Heat Flux at'. 
k ' Vorious Position, W/M**2’ . ' . '/2x , ' go. ' ) 

endi f 

9999 stop 

end 


c »«•#••*#**»•»•»»»»**••••#*#•«**•*••*•**#•***•*•*•«•*•««•*«•*****••« 

c TRIDI AGONAL MATRIX SOLVER : oWi.X(i-l) + aPi.X(i) + aEi.X(i+1) - Si 


c 

c 

c 


c 


10 

c 

c 


SUBROUTINE TDMA(X, oW.aP.oE.S, N) 

implicit double precision (o-h, o-z) 
DIMENSION X(1), oW(1),aP(l),oE(1),$(l) 
DIMENSION P(1000) .0(1000) 

( aW(1 ) - 0. , oE(N) - 0. ) 

31) : -IKJfclii 


DO 10 i»2,N 
DENUM - oP( 1 


Q 

CONTINUE 


•P( 1-1) 

( S(i)-oW(i)*Q(i-1) ) /DENUM 


P£j| « — oE( J ) /DENUM 


X(N) - Q(N) 

DO 20 i=N-1 .1 , -1 

X(i) - P( i )*X( i+1 )+Q( i ) 



